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M321 RI Broadcast Notes 
RADIO PROGRAMME 1 k 


INTRODUCING THE COURSE 


You are advised to have read Section 2 of W as far as equation (2.13) on W12 and the 
associated notes in Unit.1 thoroughly before listening to this programme. Have W on 
hand during the broadcast as extensive reference is made to it, The figures below may 
also be useful. 


Ralph Smith introduces the course, laying emphasis on the fact that the student is 
not expected to attempt every self-assessment question in each unit. Geoff Moss then 
discusses d'Alembert's solution of the wave equation as applied to a finite string and 
describes the behaviour expected of the string when it is initially at rest and has the 
shape of a pulse function (a pulse function about some point x, is a function which 
has image zero except in some small neighbourhood of x). Ralph Smith reviews the 
means by which d'Alembert's solution can be extended to take account of fixed 
boundary conditions. This procedure is applied by Geoff Moss to explain how the 
moving wave pulses generated by the original pulse function bounce off the fixed 
ends of the string. The subsequent motion of the pulses is then outlined. The pro- 
gramme is concluded by Ralph Smith, who invites students to comment on the 
course's broadcasts. 


During the programme reference is made to equations (2.1), (2.5), (2.6), (2.7), (2.8), 
(2.9), (2.11) and (2.13) of W. 


The moving points referred to in the programme are shown below: 


Moving point Moving point 

for which for which 

X + ct = constant X — ct = constant 

4——————5 -——> 
speed c speed ¢ 


Fig. I 
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If f(x) is a pulse function of height H about the point x = x, it will be of the form 
shown in Fig. 2(a); the graph of u(x,t) against x for fixed t > 0 in d'Alembert's 
solution with f a pulse function and g — 0 will be as in Fig. 2(b). 


u(x,0) = f(x), ^ 


a pulse function 
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Å pulse is approaching the end x = 0 of the string where a fixed boundary condition 
holds (u(0, t) = Ofor t > 0). Just before the leftward bound pulse hits the boundary the 
situation is as in Fig. 3(a). Just afterwards a rightward bound pulse of the opposite 
sign has emerged. This is shown in Fig. 3(b). 


Slx + et) 


fix- et) 


A 


A 


t = (x, — yc 


Here f(x + ct) = ip E ja = Y 


= f(x). 
Fig. 3(a) 
x 
t= (x, + OÓyc 
Heer x +ô 
ere f(x — ct) = f| å - e4——— 
m 
= f(-x) 
= —f(xi) by odd 
extension about x = 0. 
Fig. 3(b) 
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RADIO PROGRAMME 2 Broadcast Notes 


LINE AND SURFACE INTEGRALS 


Ralph Smith introduces the concept of a line integral and after giving a definition 
looks at a specific example. Peter Smith gives the mathematical description of a 
surface S, describes how an element of area AS can be evaluated and works through a 
surface integral problem. Reference is made to the following notes in the programme. 


[og 


Fig. 1 
Work done along P;P;41 œ f(P,)- Ar; 
«dr= li j)- Ar; 1 
[: dr = im, FAP): Ar () 
b 
[ra - I ficto) a Q) 
JC a 


Example 
Evaluate I f dr when f(x, y, z) = (y, x, 1) and C is given by 
c 


r(t) = (cost,sint,t) for0 <t < 4m. (3) 


t=4n 
5 
Cc 
1= 27 
1=0 Fig. 2 


Solution 
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f(r(t)) = (sin t, cos t, 1) 


é = (—sin t, cos t, 1) (4) 
f. T = —sin?t + cos?t + I = cos2t + 1 (5) 
4n i 4r 
f (cos 2t + 1) dt = [s + i = 47 (6) 
0 Dod o 
Exercise 
enl f- dr when f(r(t)) = (sin t, cos t, 1) and C is given by 
t)=(1,0,2) forO <t € 47. (7) 
y if 
je Fig. 3 
(a cos u sin v, a sin u sin v, a cos v) (8) 
v 
in 
0 2n u 
Fig.4 
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The integral of a scalar field ¢ over a surface S is written as 


ff ods. 


u = constant 


v = constant 


AS = area of PP,P,P; œ |PP,| PP; sin 0 
—PLI£E 


= IPP, x Pr 


= |r, — x (2 — r)| 


ôr ôr 
E A 


ôu dv 
The integral of the scalar field $ over S is 


jr 


Example 


du dv. 


or T or lu do 
Qu — Qv Mr 


Evaluate ff $ dS where @ = x? and S is the surface of the sphere x? + y? + 2? = a?. 
s 


Solution 
The surface can be represented parametrically as 
r = (acosu sin v, a sin u sin v,a cos v) O=Tu<2x0<v<n 
On the surface, 
$ = a'!cos?usin?v. 
Now 
or 


— = (—a sin u sin v, a cos u sin v, 0), 


Qu 


or å i 
— = (a cos u cos v, a sin u cosv, —a sin v), 


ov 
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so 
i i k 
êr ar "i è n 
—x—=|-asinusinv acosusinv 0 
Qu ov i : 
acosucosy asinucosv -—asinu 
= a^ [—cos u sin?vi — sin u sin?vj — sin v cos vk], 
and 
ôr ôr : å 
— x —| = a?(sintv + sin?v cos?v)? 
Qu Qv 


— a?sin v. 


We can now evaluate the integral: 


fies- [Lez 


or 
—x 
Bu dp du dv 
m pån 
= f Í (a?cos?u sin?v) (a?sin v) du dv 
0 4 


2 


n " 
«[ cost du | sin?v dv. 


tl 


0 0 
Now 
2n 27 
f cos?u du = f X1 + cos 2u) du = 7, 
0 0 
z ven 
sin?v dv = -f (1 — cos?v)d(cos v) 
0 v=0 
= — [cos v — $cos?v]yzz 
= 214 
=2—3=% 
so that 


Ii $ dS = $na*. 
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TELEVISION PROGRAMME 1 Broadcast Notes 


BOUNDARY CONDITIONS AND 
PROPERLY POSED PROBLEMS 


TV1 may be considered revision of certain aspects of the first three units of the course. 
In the programme we examine 


(a) the way in which boundary and initial conditions determine a unique solution 
of the partial differential equation modelling a physical situation, 


(b) what happens when there are too few boundary conditions, 
(c) what happens when there are too many boundary conditions, 


(d) the importance of continuity of a solution with respect to the data from which it 
is derived. 


The ideas are tackled through examples in each case. Ralph Smith introduces the 
programme and outlines its contents. ` 


Determining a Unique Solution 
Daniel Lunn discusses various physical phenomena with wave characteristics, all 
satisfying the wave equation 
2 
Uu — C Uyy = 0. 
[A subscript notation for partial derivatives is used throughout the programme. Thus 


ĝu Ou tc] 
AD. = Lletc. 

år 

He displays several solutions of this equation and explains that the correct solution 
corresponding to any given physical situation is deduced by using the associated 
boundary and initial conditions. 


y= 


This is illustrated by the case of a string fixed at both ends. The appropriate boundary 
conditions are 


u(0, t) = u(l, t) = 0. 


The string could be set in motion by being plucked (drawn aside from its equilibrium 
position and released from rest) or by receiving an initial impact while in its equilibrium 
state. The initial conditions corresponding to these two possibilities would have the 
forms 


u(x, 0) = f (x) Oxxcxl 


ĝ 

Pi 0) = 4(x,0) = g(x) Oxxxl 

respectively. In general both types of initial condition are required. The four con- 
ditions imposed on the string are then sufficient to determine a unique solution to the 


wave equation. 
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Too Few Boundary Conditions 


To show that this approach does not always work, Geoff Moss considers the physical 
situation of a long narrow ice rink of width I, insulated along one long side and 
cooled to a temperature of — 10°C on the other by a pipe carrying a refrigerating 
substance. This can be modelled mathematically by an infinite strip in the xy-plane; 
we suppose that after any initial fluctuations the steady-state heat conduction 
equation for two dimensions will hold. This is just Laplace’s equation 


Ue + up = 0. (1) 


[This is not strictly speaking the correct equation, since unless the surface of the rink 
is insulated, or equivalently the atmosphere is at a constant temperature of — 10°C, 
there will be a transfer of heat at the interface between the ice and the air above it. 
A suitable equation to describe this situation would be Uxx + Uy, = h(u — uo), where 
Uo is the atmospheric temperature and his a constant. This is difficult to solve however, 
so we imagine a simplified physical situation in which no heat transfer takes place 
at the surface of the rink.] 


The boundary condition due to the refrigerating pipe at x = l is 
u(l, y) = —10. 2) 


The boundary at x = 0 is assumed insulated. There cannot consequently be any 
flow of heat across this boundary. Since heat flow is proportional to temperature 
gradient, the component of temperature gradient perpendicular to x = 0 must be 
zero. That is 


u.(0, y) = 0. (3) 


The equation (1) together with the boundary conditions (2) and (3) has the simple 
solution 


u(x, y) = —10 a 
but we also find the further solutions 
u(x, y) = —10 + 0677/20 595 z) 5 
( — 32yj21), 3nx 
u(x, y) = —10 — 3032/2050 al (6) 


[N.B. There is an error in the TV display of equation (6), where a y appears in place 
of x on the right-hand side.] 


There are many other solutions as well. In order to find a unique one we must impose 
further conditions. We see in equation (5) that u ~œ as y ~ —co. In equation (6) 
u is in general unbounded as y ~ — co. [It is stated in the programme that u ~~ — oo 
as y ~ —oo in equation (6). In fact this is only true for 0 < x < 1/3, and u ~ co for 
3 < x < Lu = —10 for x = 1/3.] In equation (4) however u is constant and hence 
bounded as y ~ +00. An infinite temperature has no physical meaning, so we 
impose the extra “boundary condition”: 


u remains finite as y ~ +00. 


We find then u(x, y) 2 —10 to be the unique solution of the problem. The non- 
uniqueness arose from specifying too few boundary conditions. 
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Too Many Boundary Conditions 


Daniel Lunn next investigates whether il is possible to have too many conditions. He 
tries to solve a Dirichlet problem for the partial differential equation u,, = 0, i.e. he 
seeks a solution in the region R subject to its values being known at all points of the 
boundary C. The boundary values for u are chosen in the most general possible form, 
subject to the condition that u should be continuous along the boundary. R in this 
case is a square of side I. To make u continuous along C we require the functions 
fis fas 81» 82 to satisfy compatibility conditions (see picture). 


ug 
COMPATIBILITY 


The general solution of u,, = 0 is 


Ei 


x, y) = p(x) + gly), 


where p and q are arbitrary functions to be determined by the boundary conditions. 
Putting x = 0 gives 


PO) + 40) = fi) (7) 
and putting y = 0, 
P(x) + q(0) = 216). (8) 


We can then find p and q from these equations, but there remains one arbitrary 
constant in these functions (if p(x) and q(y) satisfy the equations then p(x) + K and 
q(y) — K will also do so, K being any constant). Choosing q(0) = 0 we find from (8) 
that 


P(x) = g(x) 

and then from (7) that 
ay) = fi) — 2100). 

Hence we obtain the unique solution to the problem. 
u(x, y) = gi) + fil») — 210). 


[Since u(x, y) = p(x) + q(y), the arbitrary constant K mentioned above cancels out 
in the final solution. It is not therefore necessary to choose q(0) — 0; we could instead 
just add equations (7) and (8) and substitute for p(0) + q(0) using (8) with x = 0.] 


We have apparently obtained a unique solution using only two of the four boundary 
conditions. Since the remaining two conditions express u in terms of arbitrary func- 
tions, the "solution" will fail in general to satisfy them. The problem is overdetermined 
and it is not in fact meaningful to talk of a solution. Dirichlet’s problem for the 
equation u,, = 0 has no solution. 


13 
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Continuity With Respect to the Data 


So far the need to investigate both the existence and uniqueness of solutions has been 
demonstrated. It remains to show that a useful solution should be “continuous with 
respect to the data”. Dominic Jordan defines this last phrase as meaning that if the 
initial or boundary conditions change by small amounts, then the corresponding 
solution experiences changes of the same order of magnitude. He describes it more 
loosely as “small causes produce small effects” and refers to Weinberger for a precise 
statement. 


The principle is illustrated by looking at a method for evaluating the temperature 
distribution in a fire. A poker is pushed into the fire and left to settle down. When it 
is pulled out some time later it will contain an approximate record of temperature 


distribution at the instant of withdrawal. This record will rapidly deteriorate; heat 
loss could theoretically be prevented by suitable insulation, but heat would still 
transfer from hotter to colder regions. Even so an uneven distribution would remain, 
and it might be possible to deduce the state of the fire from it. The first thing to do is to 
find out whether the problem is well posed. 


At t = 0 the temperature u(x, 0) is given by a function g(x). We suppose this to be the 
unknown temperature of the fire. At a later time the distribution is given by f (x), 
which is measurable. In the meantime we assume the poker to be completely insulated, 
so that u, — 0 at the ends and the heat conduction cquation is obeyed. We seek g(x) 


14 
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given f (x) and it looks like an initial-boundary value problem turned upside down. 
It soon becomes apparent that small changes in f will signify large changes in g; 
one way of showing this is to imagine two different fires, whose pokers show even 
temperature distributions save for a small narrow bump in one case. We infer the 
temperature distributions in the fires to be fairly even, except for a hot spot with a 


very high temperature in the fire whose poker displayed the small irregularity. It 
seems that the problem is not well posed, and that the solution will not be con- 
tinuous with respect to the data. This is significant even in the absence of hot spots 
in the fire, since the data function f (x) would in practice consist of measurements 
subject to small random errors. The computed data would be only an approximation 


to the actual temperature distribution at the time of measurement, and a computer 
programmed with a scheme to obtain g(x) from f(x) would determine a closely 
spaced series of spurious peaks. 


Dominic Jordan concludes by inviting you to turn the “upside-down initial value 
problem" the right way up by using the variable —t in place of t, comparing your 
findings with statements on W61. 


Geoff Moss and Daniel Lunn summarise the programme. We have emphasised the 
role of boundary conditions in partial differential equations and the need to know 
that a problem is well posed. Only ifa solution exists uniquely and is continuous with 
respect to the data do we have a properly posed problem, and only then can we 
model a given physical situation with any validity mathematically. 


15 


M321 R3 d Broadcast Notes 
RADIO PROGRAMME 3 


FINITE-DIFFERENCE METHODS 


The programme is linked with Unit 5, the first numerical unit of the course. After an 
introduction by Peter Thomas, Leslie Fox proceeds to demonstrate that the finite- 
difference methods used to solve partial differential equations are an extension of the 
methods for ordinary differential equations described in Unit M201 21, Numerical 
Solution of Differential Equations. He illustrates this by discretizing the onc- 
dimensional heat equation ĝu/ĝt = 0"u/öx" first in the x-direction, then in the 
t-direction. Peter Thomas discusses the variations in accuracy and efficiency between 
different finite-difference schemes, describing measures of global and local accuracy 
Finally the circumstances in which local errors fail to accumulate are examined. 
Reference is made to the following notes in the programme. 


The ordinary differential equation 
YO) + fe9y(x) = gO) (1) 


together with the initial condition y(xo) = « gives a unique solution. The first order 
Taylor approximation to the true solution is 


I(x) = y(xo + k) = y(xo) + ky'(xo) (2) 
where x, = xo + k. 
Substitute this into equation (1) to obtain 

PL = P(x) = pg + E) = y(Xo) + k(g(xo) — f (xo)y(x)]- (3) 
Now expand in a Taylor series about x; , repeat the substitution, and so on: 


GR = yea) = y(x, + k) = yo) + k{g(xı) — yE) 


Yer = 64) = DOG E E) = ye) + kele) — aay (4) 


where x, = xo +rk,r=1,2,.... 
y^ 
ar kn 
A 
x 


Xo 1 X; X, > Fig. 1 


Curve fitting to the points I pir = 0,1,2,. 


1 ..] gives i i 
true solution. } gives an approximation to the 


Turning now to partial differentia 
equation 


Ou Pu 
ôt ax? (5) 


equations, consider the one-dimensional heat 
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We wish to find values of u(x, t) in the region shown in Fig. 2, with specified initial and 
boundary values. 


u given 
u given 


—  Fig.2 


x 


€ > Fig.3 
Xo % X? x s 


Equation (5) is to be reduced to a system of ordinary differential equations. Discretiz- 


ing in the x-direction gives x; = xo + ih, i = 0,1,2,..., n. With 
Ow tj, — 24 + yy 
| ET ET HA 6 
E h? © 


equation (5) becomes 
Ôu; wv. — 2u; + Un, 
at h? 0 


This is an ordinary differential equation, because u; is a function of t only. It can there- 
fore be written 


du Ur —2u + ua 
"m ipm . (8) 


There is one of these equations for each value of i from I to n — 1. The system can be 
written 


du, 2 "EIS 
i E (9) 


ug and u, are specified for all t > 0 and u; is a function of t along the line x = x;. 


We can write this system of equations in matrix form as 
du(t) 


g TAO) i=1,2.,m—1 (10) 
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where A is the tridiagonal matrix 


-2 I 
1-2 1 
1 
he > 
1-2 1 
1 -2 
ust) Ug(t) 
us(t) 0 
u(t) is the vector * and b(t) the vector 73 
0 
un- yt) u(t) 


We next use the previous discussion of ordinary differential equations to discretize 
in the t-direction. [u consequently acquires a second suffix. At both stages of the dis- 
cretization a suffix ranging over a discrete set of values replaces a variable which may 
take any value in a closed interval. Thus when we discretized in the x-direction 
u(x, t) became u(t). This now becomes u;,, as we discretize in the t-direction.] 


Applying equation (4) to equation (9) gives 


Witir tH or 2 
Hire = Uy +k [erts re 


or 
k 
Nira = HL, + D (usi, 7 20, uias) (11) 


where the suffices can be identified using Fig. 4. 


i-lr ir [itr 


Fig.4 
Problem 
Show that use of the Trapezoidal Rule method given by 
Yea =), + Flyt + Vea) 


to solve equation (9) produces the Crank-Nicolson implicit finite-difference formula 
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The global error e; j of a finite-difference scheme is defined as 
enj = Vig uy (12) 


where U is the true solution of the partial differential equation concerned and u is 
the solution obtained from using the finite-difference scheme. The local error of a 
finite-difference scheme at a point is the addition to the global error caused by an 
application of the approximating formula at the point. This is illustrated by returning 
to the case of the one-dimensional heat equation 


ðU QU 
4r XE (13) 
approximated by the simple explicit scheme 
Uijer m Uj rcg — Ruj + uisi). (14) 
k 
whi => 
ere; = xg 
ij 
i= Lj ij [itd 
Fig. 5 


Assume that we know the true solution U at all points along t = jk and substitute 
the values into equation (14) to give 


uei Un + (Uii — 2075 + Vier) (15) 
The error in uf, , as an approximation to U; j+; is 
Ujj+1 m Ub jel = gc Ui; = HU;-1, = 2U;; + Uisa.) (16) 


By comparing equations (14) and (16) we see that the local error made in using the 
explicit scheme once is obtained by substituting the true solution of the partial 
differential equation into the finite-difference scheme. This quantity is one form of 
the local truncation error, TF 41: 


Tfi = Ui jer = Ui; = r(Ui,j = 2U,; * Ui) (17) 


for the simple explicit scheme (14). We may also obtain a different form of local 
truncation error, T; j, by using the finite-difference equation 


Upper dp Ujapj — 24 Hag 
ijt HjooUiedy ij i-Lj a8) 


k k? 


and defining 
C Ujer — Uig Uni 2Uij + Uii, 


= 19 
Tij k h? (19) 

We see that 
Tia, = kT; (20) 


for the schemes (14) and (18). 
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Problem 


What would be the equivalent relationship to equation (20) for the simple explicit 
scheme 


Ui TU js Meng — Quy jb i-i 


k? k 
applied to the hyperbolic equation 
ey pu] 
6 ox 
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TELEVISION PROGRAMME 2 roadcas otes 


CASE STUDY 


This programme and Radio Programme 5 which follows it examine pressure trans- 
mission in an apparatus used for assessing air flow in mines. The physical situation 
is modelled mathematically and a nonlinear partial differential equation results. 
An analytic solution is found to the linearized equation (since the nonlinear term is 
small it may be ignored to a first approximation). Numerical methods are then used 
to solve both the linearized problem and the original equation. The results are finally 
compared with data obtained from an experiment. Ralph Smith introduces the 
programme. 


The Physical System 


We move to a simulation tunnel in the British Coal Board's research institute at 
Swadlingcote (not “down a coalmine” as is suggested). Peter Thomas explains the 
main features of the apparatus. This consists of a long length of thin plastic tubing 
(initially out of focus in the right-hand half of the screen) with a manometer (which 
measures pressure differences) connected to one end. A tap at the other end of the 
tube is initially closed. After it has been opened the manometer will eventually register 
the pressure difference in the tunnel between the two ends of the tube, but there is a 
time lag before the levels of the liquid in the manometer re-approach a stationary 
position. 


We are interested in measuring the response time of this apparatus, i.e. the time which 
it takes the manometer to register some given percentage of the total pressure 
difference (90 per cent is taken in the programme). 


Mathematical Modelling 


Back in the studio Peter Smith begins the mathematical modelling process with a 
representation of the apparatus. In the first instance one might expect the wave 
equation to model adequately the motion of the gas down the tube, the assumption 
being that the gas travels at the speed of sound, 340 m s^ !. Such a hypothesis gives 
results which bear no relation to experimental evidence however; the response times 
are much larger than predicted. The only likely explanation is that viscosity plays 
a part due to the small size of the tube. 


Three equations are used to describe the motion of the gas. The first of these gives the 
total flux or flow rate q in terms of the coefficient of viscosity of the gas y, the radius 
of the tube a and the pressure gradient p, (as in Television Programme 1 a subscript 


21 
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notation is used to denote partial derivatives). This equation, which was derived for 
the case of steady flow in Unit 3, is 
4 


EEE 1 
q &p ^ (0) 


[It is shown in the solution to SAQ 14 on page 39 of Unit 3 that the flux Q through 
a circular pipe of radius a is given by åkna*; k here is a positive constant equal to 
p/pv where p is the pressure gradient, p the density of the fluid and v its kinematic 
viscosity (see page 20, lines —6, —5 of Unit 3). In the present restatement of the 
equation Q has been replaced by q, p describes the gas pressure (instead of the pressure 
gradient) so that the pressure gradient has magnitude |p,|, and the kinematic viscosity 
v is replaced by the (constant) coefficient of viscosity 7, where y = pv. The minus 
sign is necessitated by the fact that there will be a flow of fluid in the positive x- 
direction only if the pressure decreases with x.] 


D, zg (04). j= $ "(Constant) 


P, =D (pp).D = m za Ep 


The second equation relates the density p and the flow rate g: 
1 23 
Pr = OM (2) 


[This is the equation of conservation of mass and not the equation of motion as is 
erroneously stated in the programme; equation (1) is the equation of motion.] 


Equation (2) is derived by considering a small section of the tube of length óx. The 
mass of this section of gas is ra?óxp(x, i) so that the rate of increase of mass with 
respect to time is za? p,óx. We now obtain a second expression for this rate of increase. 
Mass flows into the left-hand end of the small section at a rate p(x, t)q(x, f) (remember 
that q is the integral of fluid velocity over the cross-sectional area of the pipe and 
consequently expresses the volume flow rate) and flows out the right-hand end at a rate 


p(x + dx, t)g(x + dx, t) = p(x, g(x, + S(p(x, t)q(x, t). 
The net rate of influx of mass is therefore 
pq — (pg + ó(pq)) = —ó(pq), 


and we may use Taylor's theorem to approximate this by —(pq),óx. Equating the 
two expressions obtained for the rate of increase of mass and cancelling óx gives the 
required result. 


Thirdly we assume that the temperature of the gas remains constant throughout the 
motion. Boyle's Law then states that the pressure of the gas is proportional to its 
density. We therefore have 


P Po 

ED G) 

P Po 
where Po, Po are the initial values of the pressure and density respectively (they are 
also the values of these quantities in the vicinity of the manometer (but outside the 
tube) for all times t > 0). 


Elimination of p and q between equations (1), (2) and (3) gives a nonlinear partial 
differential equation with p as its dependent variable, 
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Pi = D(pp,), (4) 


where D = a?/8y is a constant (not to be confused with any form of differential opera- 
tor). 


The transformation 
P — Po 
P= 
Apo 


where Ap, is the pressure difference in the tunnel between the two ends of the tube, 
is now used to simplify the problem. e is written for the ratio Apo/po. Equation (4) 
becomes 


P, = «(1 + 6P)P,], (6) 


where x = a*po/8y is a constant. 


(5) 


Since p = po at t = 0 we see from (5) that 
P(x, 0) = 0 D<x<I (7 


where I is the length of the tube. After the tap is opened at x = I the pressure there is 
maintained at py + Ap, so that 


PL)-1 t>0. (8) 


At the manometer we assume that the flow of gas is zero, i.e. q(0, t) = 0. From equation 
(1) this means that p,(0, t) = 0, and so from (5) 


P(0,:)—-0 t>0. (9) 


P=0,/=0,allx - 
P=1,00,x=l — 
P=0,0>0, x=0 


Equations (6) to (9) give a problem composed of a nonlinear partial differential 
equation together with initial and boundary conditions. In practice ¢ is found to be 
small so that it may be neglected to a first approximation. Equation (6) is then 
replaced by 


P, = KP, (10) 
which is the heat conduction equation. 


Equations (7) to (10) give a standard boundary value problem. The viewer is referred 
to the broadcast notes for its solution. [In fact the solution to this problem is derived 
in Radio Programme 5 and can be found in the corresponding broadcast notes.] 
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What Does the Model Predict? 
At the manometer we find 


E 209 (=o? iia |. 
P00, =1 2E reno (n+4)n E (11) 


P(0, t) approaches the value I as t becomes large. Suppose t = tg when P(0, t) reaches 
some fixed percentage of its final limiting value, say 90 per cent for definiteness. 
Then P(0, tg) = 0.9, and we can find a constant c such that 


0921— : Y ET exp[—(n + 3?z?q. (12) 


zon + 


be 


It is clear that c is independent of the dimensions of the tube and of the response time. 
Putting t = tr in equation (11) and comparing with (12) we see that ktg/” = c. 
Using the expression for x below equation (6) it follows that the ratio a?ts/P is 
constant as the length / and the radius a are permitted to vary. In other words the 
response time tg is proportional to the square of the length of the tube and inversely 
proportional to the square of its radius. 


The Numerical Approach 


Peter Thomas looks at the linearized problem given by equations (7) to (10) with a 
view to solving it numerically. A piece of computer-animated film shows the dis- 
cretization of the rectangular region in which a solution is required and of the initial 
and boundary conditions. [Two points here: we cannot “take the response time to 
be T” as stated in the programme, since it is the response time which we wish to find 
by solving the problem. T should therefore be chosen to be considerably larger than 
any expected response time. The film successively assigns the values 1 and 0 to Py o. 
In practice Py,o is put equal to 0.5 rather than either of these values (see the notes to 
Library Program $RTT 321). We then have Py; = 1 only for j = 1,2,...,M, and 
not “from 0 up to M” as stated. P incidentally is being used here as a symbol for both 
the analytic and finite-difference solutions for the pressure difference.] 


TM| 


h and k are taken to be the (constant) mesh sizes in the x- and t-directions respectively, 
so that Nh — I, Mk — T. The derivative boundary condition at x = 0 is dealt with 
by using the central difference formula 


(Po; = (Pi; = P_, )/2h (13) 


and extending the mesh accordingly. The discretized initial and boundary conditions 
are 


Pi9 =0 i= -1,0,1,...,N—1 
Pyg = 0:5 
Pyj=1; Paajy=Pi; j=1,2,...,M. (14) 


Equation (10) is then replaced by the finite-difference scheme 


Pipi ~ Pig Pay yt Pip 
k n n (15) 
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Putting the mesh ratio k/h? equal to r this can be rearranged as 
Bgm Pj Kr(P. c; — 255 P4) (16) 


PE P. 


Pj uS P; HC Y(Pj.,,-2P,; +P; 


Os rss 


which is an explicit scheme. It is only stable when r satisfies the inequality 0 < xr <4, 
and this condition is somewhat prohibitive in terms of practical computing as is 
shown with a numerical example. The stability condition is k < h?/2x. Putting 
« = 100,000 (a typical value), | = 100m and N = 5, we have h = 1/5 = 20 m, giving 
the restriction k < 0.002 s. Even if we take the equality here a value of M = 1,250 
is needed in order to reach T = 2.5 s. Though this is a large value for M the situation 
becomes progressively worse if we seek a more accurate scheme. 


If ten times as many grid points are taken on each time row we require M to be 
125,000 (not “13 million” as stated in the programme, but still too much for a computer 
to handle). 


This stability condition is avoided by using an implicit scheme to replace equation (10) 
(for example the Crank-Nicolson scheme 


—KrP iy ig, + (2 + 2Kr)P; 541 = Pi pay = KrP; Lj 
17 
+ (2 — 2xr)P; j + KrPi4 sg a7) 
is à satisfactory replacement, being stable for all values of 7). In the nonlinear case 
the situation is very similar. A typical value of 0.01 is taken for c. The initial and 
boundary conditions remain the same. An implicit scheme is again employed to avoid 
stability restrictions on the mesh size (although the viewer is referred to the broadcast 
notes for details of the implicit method in the nonlinear case this again is dealt with 
in Radio Programme 5). 


Comparison With Experiment 


We return to Swadlingcote to conduct an experimental test of our results. The pressure 
difference, created by a large fan, is known in advance and the manometer is marked 
at a point representing 90 per cent of the final reading. The time which it takes the 
liquid in the manometer to reach this level after the tap has been opened is recorded 
on a stopwatch. 
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Finally a set of such data for various different lengths of tube is compared with values 
obtained from the analytic and numerical solutions. The line on the left of the picture 
joins the experimental results while that on the right represents all three methods used 
in the programme; it is impossible to differentiate between the methods on this scale. 
[You may well wonder why these graphs are straight lines when it was shown earlier 
that at least in the analytic case the response time is proportional to the square of the 
length. The explanation is that the scales are not linear but logarithmic (this is the 
significance of the unequal mesh spacing on the graph paper). We can alternatively 
consider the graph as plotting log tr against log! on linear scales. Since tg cc P, 
log tg = 2log I + constant, giving a straight line] Since each method examined 
produced a result close to that of the experiment it is concluded that the linearization 
of the problem was a reasonable approximation to make. The viewer is referred to 
Radio Programme 5 for further discussion of the case study. 
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RADIO PROGRAMME 5 Broadcast Notes 


CONTINUING THE CASE STUDY 


This programme follows on directly from Television Programme 2. After an intro- 
duction by Peter Thomas, the physical aspects of the case study and its mathematical 
modelling are described again by Peter Smith, who subsequently solves the resulting 
(linearized) initial-boundary value problem using separation of variables. Peter 
Thomas then discusses the numerical solution of the nonlinear case, introducing a 
suitable implicit finite-difference scheme and describing how problems which arise 
in attempting to implement the scheme as a computer program can be overcome. 


The paper by Jones and Jordan which is referred to in the programme is reproduced 
on the pages following these notes. This is followed by a listing of the library program 
$RTT 321, notes as to its use and a set of exercises. Reference is made to the following 
notes in the programme. 


Figure I is a diagramatic representation of the apparatus. 


Fig. I 


The volume flow rate q in the pipe is given by 
på EE. q) 


where a is the radius of the tube, 1 the coefficient of viscosity and Op/Ox the pressure 
gradient. 


We assume that the gas satisfies two further equations, 


op 1 apa) " 
Ot — na? ox 

and 
P_ Po, (3) 
P Po 


which are deduced from the law of conservation of mass and Boyle’s Law respectively 
(see the TV2 broadcast notes for more details about this and the following). The initial 
and boundary conditions for the pressure are 


p(x,0) = po O<x<l, 
Pt) —po-c Apo t>0, 


Poy 9 t0. 
Ox 
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After eliminating p and q between equations (1), (2) and (3), and replacing the pressure 
p by the pressure difference P = (p — po)/Apo, we arrive at the nonlinear partial 
differential equation 


op a oP (4) 
9^ “Ze + eP) z] 
where £ = Apo/po and K = a?po/85. Neglecting the nonlinear term, 
oP 0?p (5) 
a ae 
4 


Fig.2 
For convenience P is replaced by U — P — 1, so that, 
U(Lt) = 0 t>0, 


oy 0 t>0, 
ôx 


U(x 0 = —1 Oxxx«l 
Put 

U(x, t) = F()G(t), (6) 
so that F and G satisfy 

FG s 

F xG ^ (7) 
for some constant c. F satisfies 

ær >, 

dé +cF=0 (8) 


and the general solution of this is 
F(x) = A cos cx + B sin cx. (9) 


The boundary condition at x = 0 implies that B = 0, and the condition atx=| 
implies that cos cl = 0. c can therefore take any of the values n/2l, 3n/2l, Sn/2l,... 
or in general (2n + 1)/21 for integers n. Substituting this into the equation for G, i 
dG Qn + 1)?x? 
[MET ee (10) 


This has the general solution 


» (2n + 1yz 
G= De | «PET, (11) 
where D, is a constant. 
Adding all the separable solutions together we get 
= —x(2n + 1?z?t (2n + 1) 
Uj) = Y D,exp| e zx 
) x p | 7P eo M (12) 
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and the D, are chosen so that the cosine series obtained by puttingt = Oin equation (12) 
is equal to —1 on the interval [0, I]. The required series is 


e (c1) (n + l)æx 


adm ENE 
olan + a S 3) 
so that the full solution is seen to be 
2 (=1F —k(2n + 1?z?t (2n + lax 
P(x,th=1-4 d 
) Xm. eo? ra COS 5i (14) 


We find the pressure difference at the manometer by putting x — 0 in this equation. 
Note that the time t always appears in a factor xt/P”, so that this factor must always 
be the same value (for different lengths and radii of pipe) when P(O, t) achieves any 
Biven value. The response time is therefore proportional to the square of the length of 
the tube and inversely proportional to the square of its radius. 


For equation (4) there exists an implicit numerical scheme which is both stable and 
accurate but gives rise to systems of linear algebraic equations, namely 


Pos — Piya = Sefo (Phi jei Pipe) + (Pisiy PL) 
od (Pi+1,j-1 = Pi j-1)} = a7 (P; = Pia jaa) (15) 
+ (Pig — Pina) + 0,24 — Pi-ay-)) 
where 
at = 1 + de(Pi4.; + Pih 
a” = 1 +3e(P; PC) 


andr = k/h? (P is being used here as a symbol for both the analytic and finite-difference 
solutions for the pressure difference). 


The initial and boundary conditions for the problem in terms of P are 
P(x,0 =0  O<xx<1, 
Pt =1 :»0 (16) 


The discontinuity in the data at (1,0) may affect the numerical solution, but this can 
be minimized by choosing a very small time step initially. Since equation (15) repre- 
sents a three-level scheme we need to know values along / = 0 and j = 1 before we 
can start using it. Values for j — 1 are calculated from the known values on j = 0 by 
using a simple two-level scheme. We can construct a suitable explicit scheme as 
follows. 


Writing equation (4) in the form 


ôP ap? 0?P 
== z P)—; 17 
à [ (FI AS aE (17) 
and replacing the first and second space derivatives by central differences, 
oP ~ Pai Pinta, (18) 
Ox]ij 2h 
P " Pia — aP i,j + Prati, (19) 
x] i; h? 
we arrive at the scheme 
Kr 
Pipa = Pij + 4 [E(Pi+1,5 — PL) 
(20) 


+ A(L + EP MPi+1,5 — 2Pig Pici gl. 
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Since the object of using an implicit scheme is to avoid the restrictions on step ue 
which accompany explicit schemes we use the explicit scheme several times in order 
to provide values along the first time level of the implicit scheme. 


[In fact we use the explicit scheme to provide values along the first and second time 
levels of the implicit scheme so that the latter escapes the immediate impact of the 
discontinuity in the initial and boundary data at (1, 0).] 


Explicit 


Implicit 


The step length k, of the implicit scheme must be an exact multiple of the step length 
kg of the explicit scheme. The implicit scheme gives rise to a tridiagonal system of 
equations which can be solved by the recurrence method of Unit 5. 


The following paper is reprinted from the British Journal of Applied Physics, 13, 
420—423, Aug. 1962 by permission of the Institute of Physics. 
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Reprinted from the British Journal of Applied Physics, VOL. 13, pp. 420-423, AuGust 1962 


Time lags in the transmission of pressure disturbances 
along long lengths of small bore tubing 


C. JONES and D. W. JORDAN 
Mining Research Establishment, National Coal Board, Isleworth, Middlesex 
MS. received 11th April 1962, in revised form 23rd May 1962 


When it becomes necessary to transmit small air pressure differences over long distances 
the effect on the response time of the small bore tubing used to connect the source of 


pressure to a manometer should be known. 


This article develops a theoretical treatment 


for evaluating the time and shows it to vary directly as the square of the length and inversely 
as the square of the internal diameter of the tubing. The theory is confirmed as regards 
variation with length by laboratory experiments with Neoprene tubing having a bore of 
is in. The experiments also include an examination of the effect of varying the magnitude 
of the pressure, the position of the manometer and the coiling of the tube. 


1. Introduction 


The transmission of differential air pressure to a distant 
position such as occurs in the remote control of fan speed 
or the recording at a distance of the pressure drop across 
an orifice plate are instances where long lengths of small 
bore tubing are used. The question arises as to how the 
tubing affects the response time of the manometer when 
lengths of tubing occasionally as great as one mile are used. 
Of special interest are the effects of the length and diameter 
of the tube on the time. This paper describes investigations, 
both theorctical and experimental, into this question for the 
case when small bore tubing is connected to a liquid-in-glass 
manometer whose volume change has been neglected. 


2. "Theoretical treatment for the propagation of pressure in a 
Jong small bore tube 


Jt is assumed that the flow of air into the tube must, of 
course, be initiated by a wave motion, but as this is damped 
out relatively quickly compared to the times observed the 
Characteristics of the air motion are not governed by the 
propagation of sound waves, shock waves or the like but 
by laminar viscous flow of air. It is shown later that the 
results obtained are consistent with this assumption. 

The well-known Poiseuille expression for laminar flow is 


om 


mat Op 
89 ox 
where g is the volume flow rate in direction x, a the radius 
of the tube, 7 the dynamic viscosity of the air and p is the 
pressure, assumed in deriving this formula to be uniform 
over the cross section. 

A second relation connecting pressure and flow rate is 
obtained as follows. i 

Consider a short portion of the tube between x and x + dx. 
The mass flow rate at x is pg, where p is the density at x; 
and at x + 8x is pq + 5 (pg). The net influx of mass into 
the volume za?8x per unit time is therefore —{0(pq)/dx}8x, 
and the rate of influx of mass per unit volume, which is 
equal to the rate of change of density, is —(1/2ra?){0(pq)/dx}. 
Thus 


q=— 


39. 1209 D 
dr na Ox 


420 


With isothermal conditions Boyle's law holds: 


£f const = Pa (3) 
P Po 
where po and pg are reference values, say under the undis- 
turbed conditions. Equation (2) can then be written in 
terms of pressure and flow rate: 


3p I og) 
WT aa? Oe o 
and this is the second equation required. 
Substitution into (4) from (1) gives 
dp ydy dA >, Dr 9p 
wo (zaal Pa) Ds zx) 6) 


where D = a2/87. 

Consider that the end x = / of the tube is subjected to 
sudden change of pressure from pg to pg + Apo. The other 
end, at x = 0, has a manometer attached of negligible volume 
so that the flow is always zero there, a condition which 
interpreted in terms of p by equation (1). Then the initial 
and boundary conditions to equation (5) are 


p — pogatt =0, forallx 

p pod Apo fort>0, atx-! (6 
39g for?t0, atx=0 

ox 


Equation (5) is like the ordinary heat-conduction equation 
ina rod, with a variable diffusivity equal to Dp. The boundary 
conditions are analogous to a sudden change of temperature 
at one end, the other end being insulated. Equation (5) 
non-linear, but the non-linearity can be removed by demon- 
strating that its effect is negligible. Putting 


-PZ Po). Apo _ P 0 
Apo Po 
then 
oP 3 dP 
— =K = 8 
3 «lu +EP) x) (8) 
with 
P-0atr-0, forallx 
P=1fort>0, atx=/ (9) 
3P o for i2 0, atx =0 
Ox 
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where 
apo 


k = Dp, = Em 


(10) 
Ap, will be of the order of a few inches of water gauge at 
most, so that e, being the ratio of this quantity to atmospheric 
pressure, will be of order 175, and it is intuitively clear that 
P will not exceed unity. The factor 1 + eP in (8) is therefore 
effectively 1, and this makes the equation linear, its form being 


dP_ WP 

va 

The solution to this problem is given by Carslaw and 
Jaeger (1948, p. 87): 


(D 


P(x, 1) = ; = (—1)" cos FÆSE — exp (~ rpp NY Br 


2 "AS T 
It can be verified that i £ (— 1)" (cos B, x)/8, is identically 
equal to I so that ase 


5 oe exp (—B,/) cos B,x. (12) 
n=0 n 


In the above B, = (2n + 1)7/2l. 
It is desired to find how long it takes for a change in 


pressure Ape at x= Í to reach x —0, At x =0 the 

expression becomes 

P0,n=1-2 8 OD ap tn + pte} 13) 
Sr ola 4 5 ds å 


The time taken for the pressure at x = 0 to attain a given 
fraction of its final value in cases where different tube 
diameters and lengths are considered is given by the relation 


kili? = constant 
or : 
length 

diameter” (14) 

Alternatively, dimensional arguments give this result directly. 

Next, consider the actual times of propagation under 
circumstances similar to the tests described later. 

For air, 7 = 1:3 x 1075 ]b ft-! sec), a = 1/128 ft with 
the Neoprene tube available, and py = 6-8 x 104 pdl ft-?, 
so that x = 4 x 10*. Likely values of l lie between about 
300 and 10000 ft, and the experiments that are described 
later in the paper give values of 90% response time between 
about 2sec and 2 x 10*sec. Calling this time Ty:9, the 
value of «Tp.9//? has nearly the same value for all the 
experiments, as expected, the mean being about 0-90. For 
this value the first term of equation (13), corresponding to 
n = 0, has most significance, giving very nearly 


— 4 
P= EU a 1 == exp {— brer) (15) 
for values of «¢//? greater than about 0-1. If E represents 
at any instant of time that fraction of the total expected 
increment Ap, which still has to appear on the manometer, 
then E = 1 — P = 1-27 exp(— 2-46kt[P). 

It has been mentioned above that in these experiments 
xij? = 0:90. This corresponds to a theoretical value of 
E of 0:15 and implies that, if the theory were correct, the 
times were recorded when the manometer reached 85% of 
its final value instead of 90%, which can be considered a 
very good agreement. 
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That the flow rates obtained above are consistent with the 
model used, which assumed laminar viscous flow, can be 
shown as follows. The Reynolds number R, which deter- 
mines the type of flow, is given by R = 2avp/n where p 
the density and v the speed of flow, equal to g/za^. Also 
q is given by equation (1), so that 

= ap [òp] 

~ dn? [ox 
Now [dp/dx| has its maximum value at x = I, which is the 
open end of the tube. This can be seen by considering that 
dp/d1 is obviously of the same sign everywhere, and so there- 
fore is d?p/dx?, from equation (11). But the value of [ap/dx| 
is zero at x = 0, so its value increases steadily up to x =|. 
To find its value approximately at x — /, consider the case 


i (16) 


‘of a semi-infinite pipe extending from x = / to x = — oo, 


since near / = 0, when the gradients are steepest, this case 
does not differ appreciably from the case of a finite pipe. 
The solution for such a pipe is (Carslaw and Jaeger 1948, 


p. 43) 
x—i 
p = Ápy| 1 — erf Xy 


3p Am 
Ox (me)? 


andat x =] 
ay) 


The Reynolds number decreases to a value of about 10? 
appropriate to laminar flow, in a space of time given by 


=P AP _ igs 
T 47? [CHE = å 
ap A 1 \? 
ort «(å C i] & (2:53 x 107 9A7p, sec 


for 7% in. bore tube. Here Ap, is measured in pdl ft-2. This 


can alternatively be written 
t > (6-08 x 10-44 sec 


where W is the change in pressure in inches of water gauge. 
The flow will therefore be laminar everywhere for all but å 
negligible interval of time, since W will not generally be å 
large number. 

In this theory the volume change in the manometer has 
been neglected as is justified for a liquid-in-glass manometer. 
This could be taken into account if necessary by taking in 
place of the boundary condition at X =0 given in (9) a 
boundary condition which involved coupling with the 
equation of motion of the manometer. 


3. Laboratory measurements of response time 


As an experimental check of the law that the response 
time varies as the square of the length and to determine 
values of the time for a given diameter of tube, several 
Series of measurements were made using a water-in-glass 
manometer. One limb of the manometer was connected by 
different lengths of Neoprene tubing of in. bore, up to 
the maximum length available of 1200 ft, to a source of 
Pressure, the other limb being open to atmospheric pressure. 
This use of a manometer arises with pressure surveys along 
coal-mine Toadways. The Pressure was both applied to and 
removed from the measuring system using a two-way tap: 
one leg of which was connected to a static pressure tube 
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inserted into a ventilation duct, the other leg being open to 
atmospheric pressure. The manometer was thoroughly 
cleaned Íree of grease and dirt before and during the 
experiments. 

The first observations were carried cut with a constant 
pressure of 3 in. water gauge and the time taken for 90% 
of the change in pressure, from zero to 3 in. and from 3 in. 
to zero, was noted for various lengths of tube. Table 1 
shows the results obtained, each entry being the average of 
five observations. Whether the pressure is applied or removed 
the response times are the same. The least squares analysis 


of the results gives t = (//220-2)?:07 when I is in feet and 
t in seconds. 


Table 1. Response times for different lengths of 
tubing with constant pressure (3 in. w.g.) 


Length of tubing Time to reach 90% of 

(ft) pressure change (sec) 

Pressure Pressure 

applied Temoved 
300 2:0 2:1 
450 41 4-0 
600 T3 7-0 
750 14-2 14-2 
900 19-8 19:4 
1050 25-6 25:4 
1200 31-7 31:4 


A further short series of observations was completed by 
applying different pressures, up to 8 in. water gauge, to a 
constant length of tubing, namely 1240 ft. The mean value 
of the observed 90% response times was 38 sec, which is 
close to the value of 36 sec to be expected from the above 
least squares analysis. There was no significant change in 
the time for different values of the pressure. 

Table 2 gives the response times with one limb of the 
manometer connected to a source of pressure through 797 ft 
of the same tubing and also with the manometer moved to 
other positions along the tubing, namely 620, 400 and 0 ft, 
the rest of the length of tubing being attached to the other 
limb and at atmospheric pressure. Each tabulated figure is 
the average of 5 observations. Movement of the manometer 


along the tubing reduces the response time in accordance 
with the square law relationship, the open end of the mano- 
meter remaining at atmospheric pressure even though it has 
an increasing length of tubing attached to it. 


Table 2. Variations of manometer position along the tubing 
(pressure 0-56 in. w.g.) 


Tubing length 90% response time (sec) 
(fh) Pressure Pressure 
applied removed 
797 14:5 14-5 
620 9-0 9-0 
400 35 3-5 
0 Instantaneous Instantaneous 


Although it seems unlikely that manometers with tubes of 
unequal length would be used extensively, a difficulty arises 
on which comment should be made here. If the pressures 
on both limbs change, the 90% response time will depend 
upon the actual, but unknown, values of the pressure changes 
as well as on the lengths of the tubes. This is illustrated by 
an extreme case, when there is a large change of pressure 
on the short tube and a small change on the long tube, the 
time to read 90% of the applied pressure difference being 
governed by the length of the short tube. If the pressures 
are interchanged the time is governed by the long tube. 
However, it is generally the case that the response time is not 
greater than that appropriate to the long tube. For 
manometer with equal tube lengths this difficulty does not 
arise, since the ratio of the pressures at the manometer 
each limb at any instant is the same as the ratio of the applied 
pressures. 

Finally, the response times were measured at different 
pressures up to 8 in. water gauge using 600 ft of tubing, 
one case with the tube straight and in the other with the 
tube coiled but nevertheless free of kinks. Again there was 
no difference in the average values of the response times. 

The figure has been drawn from the data given in table 
so that the reader can easily read the response times for any 
desired length up to 10000 ft of a number of different 
diameters of tubing. 


| | "Á 
10000. bet 
8r H [- ae ref in 
6 Se LAT 
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90% response time (sec) 


Response time for different diameters of tubing. 
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4. Conclusions 

For liquid-in-glass manometers the following conclusions 
are made: 

(i) The 90% response time for long lengths of tubing 
connected in the open-end manner to a liquid-in-glass 
manometer varies as the square of the length and inversely 
as the square of the diameter of the bore of the tubing. The 
value of the 90% response time for 600 ft of 34 in. bore tube 
is 8-0 seconds. 

(ii) The response time does not depend on the value of the 
pressure or on whether the pressure is applied to or removed 
from the manometers. Neither does it matter whether the 
tubing is coiled or straight, provided it retains its natural 
cylindrical shape. 
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(iii) When a manometer has each limb connected by equal 
lengths of tubing to a source of pressure difference the 
response time is that for either length of tubing. 
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LIBRARY PROGRAM—SRTT 321 


3 Lib ary Program 


The following isa print-out of the library program $RTT 321, which is used to produce 
a numerical solution to the partial differential equation arising from the case study 
of Television Programme 2 and Radio Programme 5. The listing is followed by 
explanatory notes and a set of exercises. 


18 
2 
368 
40 
50 
69 
79 
80 
96 
109 
118 
122 
136 
148 
158 
166 
176 
189 
19Ø 
2800 
210 
226 
230 
249 
25Ø 
260 
276 
280 
230 
300 
316 
329 
336 
340 
359 
369 
378 
380 
399 
400 
410 
426 
430 
440 
450 


PRINT " MANOMETER RESPONSE TIME" 
PRINT " |  ----------2---.-2-..--.--.-- " 
PRINT 
PRINT 
PRINT 
PRINT 


U*X*INPUT*x'" 


"INPUT LENGTH OF TUBE (FTO "3 


INPUT 

PRINT 

PRINT 
INPUT 
PRINT 
PRINT 
INPUT 
PRINT 
PRINT 
INPUT 
PRINT 
PRINT 
INPUT 


L 


"INPUT BORE OF TUBE (INS) 
A 


"INPUT PARAMETER EPSILON 
E 


"NUMBER OF MESH POINTS 
N1 


"TIME STEP FOR IMPLICIT 
K 


PRINT 

READ Z»P® 

DATA .0002135,638000. 

A=A/24 

K1-2A*A*PO/ C3*Z) 

P1=3.14159 

DIM OC5ØIS»PL5Ø1-QL[5Ø]»-ALSØ1»8[5Ø]-»CC[5Ø1]8-DC5Ø1-WCL5Ø1.SC5Ø] 
PRINT ''*x*''; 

PRINT "PARAMETER KAPPA ="5K13 

N2=N1-1 

N3zN1-2 

REM 
REM 
REM 
H=L/N2 
K9=H*H/ (2*K1) 
REM 
REM 
REM 
M=INTCK/K9) +1 
K9=K/M 
REM 

REM R= 
REM 
R=K9*A1/C4KH*H) 

PRINT "TIME STEP FOR EXPLICIT METHOD ="3K9 


Hs STEP LENGTH IN X-DIRECTÍON 


M= NUM3ER OF STEPS WITH EXPLICIT SCHEME 


+25 TIMES KAPPA TIMES MESH RATIO (K/H12) 
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46Ø PRINT | " 3 sp; 
47Ø PRINT "FOR EXPLICIT METHOD ONLY TYPE E OTHERWISE [YPE ; 


480 INPUT AS 
490 PRINT - Sd 
508 PRINT "FOR COMPLETE PRINT-OUT TYPE P OTHERWISE TYPE N'3 
510 INPUT FS 

526 PRINT 

530 PRINT "**WORKING KK" 

549 PRINT 

558 IF F$z'"N" THEN 646 

569 PRINT "SOLUTIONS GIVEN AT X=ø" 

578 PRINT 


58Ø PRINT "TIME NUMERICAL ANALYTICAL ‘ERROR 
596 PRINT " SOLUTION SOLUTION" 

600 PRINT 

610 REM 

620 REM INITIAL CONDITIONS 

630 REM 


640 FOR I=i TO N1 
650 OCI1=Ø 

668 NEXT I 

67Ø  OLNIJ2.5 


680 REM 

690 REM EXPLICIT SCHEME 
7800 REM 

7180 JO 

720  J-J*1 


739 PL11-0L11*B*R«*CI*E*OCIJ2*COC2J-0C112 

740 PCN1)=1 

759 FOR I=2 TO N2 

769 PCII=COCI+11-OCI-11>*C0CI+17]-0OC1-10) 

779 PUII=OCII+RACEAPCII+HAKCL+EK*OCIIDKCOEI+L I-2KOCI3+OCI=19>> 
782 NEXT I 

798  T-K9*J 

820 FOR I=1 TO N1 

810  OLIJ-PCIJ 

820 NEXT I 

838 IF F$="N" THEN 86g 

840 gOSUB 168Ø 

850 PRINT TsPC11-S1>S1-PC1] 

860 IF P[1] >= .9 THEN 1489 

878 REM 

8890 REM HANDOVER FROM EXPLICIT TO IMPLICIT SCHEME 
890 REM 

900 IF A$="E" THEN 729 

918 IF J24 THEN 949g 

9288 IF JsM*2 THEN 989 

330 GOTO 729 

940 FOR I-1 TO N1 

950 Q@C1)=0C1) 

960 NEXT I 

972 GOTO 796 

980 FOR 1=1 TO N1 

990 OCII=QL1) 

1900 NEXT I 

1610 PRINT 

1920 PRINT "EXPLICIT SCHEME N ER" 
re E FINISHED AFTER 
1049 REM 


3 J"STEPS" 
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1858 REM IMPLICIT SCHEME 

1960 REM 

1270 J=2 

1080 REM 

19980 REM R= TWO-THIRDS KAPPA TIMES MESH- 

1199 REM å EE 
1112 Rz2*(OKI*N2«N2/ C3*LxL) 

1120 JsJ+1 


1132 A3=R*(2+E*(PC1]+PC21))> 
1140 BC1]=1+43 

11580 CC1]=43 

1169 DE117=0C11+A3*(PC21+0C21-PC1J-0C13> 
117Ø FOR I-2 TO N3 

1188 ADIZR*CI-EK.BS&CPCITAPCI-1)2) 

1199 AZ=R*C(1+E*eSKC(PCIHLI+PCII)) 

1209 A3=A1+A2 

1218  ALI3sA1 

1228 BCII=1+4A3 

1230  CtIJ-a2 

1246 DCIJEOCIIJ«AL K«CPCI-13*0CI-132 ABKCPC I1 340£ L1 2 -A3XCPC LJ40 C 
1258 NEXT I 

1260 AL=R*C14+E*.5* CP ON27+PIN31)) 

1278 A2=R*(1+Ek.5*(1+PCN2J)9) 

1280 A3=A1+A2 

12980  A[N21-A1 

1300 BLIN2)]=1+43 

1318 DCN2J-OCN2]J*A1*CPCN3J1-0CLN32J) *3kA2-A3* CPCN2 2-0 C423) 
1320 T=J*K 

1330 FOR I-1 TO N! 

1348  OLIJzPEIJ 

1358 NEXT I 

1368 30SU3 1530 

1376 IF F$="N" THEN 1406 

1388  GOSUB 1680 

1396 PRINT TsPC1),S15S1-Pce1] 

1406 IF PC11]<+9 THEN 1128 

1410 T=J*K 

1428 PRINT "IMPLICIT SCHEME FINISHED AFTER"3JS"STEPS" 
1430 PRINT 

1440 PRINT 

la50 PRINT "x**SULUTION**" 

1460 PRINT 

1478 PRINT 

14860 PRINT "RESPONSE TIME ="3T3''SECS’ 


1493 STOP 

1500 REM 

1510 REM TO SOLVE A TRIDIAGONAL SYSTEM OF EQUATIONS 
1520 REM 


1530 WC11=8C1) 

1540 S$C1]=DC1) 

1558 FOR I=2 TO N2 

1569 Z=ACI1/WCI-1] 

1578 WCI1=B8C1I-CCI-1)*Z 
1586 SCII=DCII+SCI-1I*«z 
15990 NEXT I 

1600 PtN221-SCN21/WCN2] 
16108 FOR I=N3 TO 1 STEP -1 
1620  PCIJSCSCIJ4CCIOKPUIT112/WELI] 
1630 NEXT I 

1640 RETURN 
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1656 REM g 
1660 REM TO CALCULATE THE ANALYTICAL SOLUTION TO THE LINEAR PROBL 
1670 REM 

1688 Si=ð 

1690 N=Ø 

1708 S2-(- I2 TNAEXP C GE EO KCN 50 XP LEP Le KT&T/ CLL22/ CN 5 

1710 S1=51+52 

1728 IF A8SCS8)«1.E-07 THEN 1759 

1738 N=N+1 

1740 SOTO 1700 

1750 S1=1-S51*2/P1 

17690 RETURN 

1770 END i 

Program $RTT 321 


1 Objective 


To solve numerically the problem 


oP ô ðP 

I eee == x«l 0, 
à «x teen) O<x<l t» 
P(x,0) =0 Osxxl 

Ph) =1 t» 0, 

P 0, =0 t>0. 

Ox 


In particular, to find the value of t for which P(0, t) first becomes greater than 0.9. 
This value of t is known ås the response time. The problem models a long thin tube of 
length I which transmits a differential air pressure down its length. The finite-difference 
method is used and the region is covered with a rectangular mesh. 


2 Input 

(i) length of tube (ft) program symbol L 
(ii) bore of tube (in) A 
(iii) parameter epsilon E 
(iv) number of mesh points in x-direction N1 
(v) time step for implicit method K 


3 Other Data (input via a DATA statement) 


(i) coefficient of viscosity of air (1) program symbol Z, value 1.3 x 1075 
(ii) pressure under undisturbed conditions (pọ) symbol P6, value 6.8 x 104 
4 Options 


(i) Use either an explicit or an implicit scheme (Type E or D. 


(ii) Print out response time only or get fuller print-out (Type N or P). 
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5 Intermediate Calculations and Informative Print-out 
() Evaluation and print-out of parameter x (kappa) given by 
K = a?po/8n, 
where a is the radius of the tube. 
Gi) Calculation of optimum time step for explicit scheme. 


(See Section 9, Implicit Scheme.) 


6 Initial Conditions 
If there are Ni mesh points in the x-direction then 
Pio =0 i=1,2,...,NL — 1, 


Pio = 0.5 (where P is now being used as a Symbol to describe the 
finite-difference solution for the pressure difference as well as 
the analytic solution). 


(See exercise 6.) 


7 Output 


If full print-out is required then the numerical solution of the differential equation is 
printed out at the point x = 0 at each time step together with the analytical solution 
to the linear equation (e = 0) at the same point. The difference (ERROR*) between 
these two values is also given. The final output is the response time. 

8 Explicit Scheme 


The explicit finite-difference replacement of the nonlinear equation is obtained from 
the form of the equation 


oP OP\? a 
UU cfe (22) +(+ Pas] 


from which we obtain 

Pi jet = Pit RE(Pi41,5— Pia) HA REPLY 1j—2P ig + Pina MÅ 
where 

R = xk/(4h?). 
For stability we require xr < 4 where r = k/h?. 


(See also Section 9, Implicit Scheme.) 


9 Implicit Scheme 


The implicit scheme for the equation 


OU | atus 0 (1) 
oe e zn a(U) 2-0, b>0, 


is given by 
bl jay — vin) = Srl (eager — s) Qi ti) 
sui 7 20) — € nsi 7 Mia jad) 
(un; — iau) + (jr 7 i-i-a) 


where 


Ui jg + Ui = Us; t Uig 
E os o a md oe 


and r is the mesh ratio k/h?. The scheme is a three-level scheme since it involves values 
of u on the three time levels j + 1,j andj — 1. 
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For the problem 


10P a zr) 

BEE P)— O<xx<1l, t>0 

x t hrag 

P(x,0) =0 O<x<i, Q 
P) =1 t>0, 

P 0) =0 t>0, 

Ox 


we see that in the notation of equation (1), a(P) = 1 + £P, b = lk, and so 
at =1+ iP + Pi» 
a^ = 1 + 4e(P,; + Pii 


It is more efficient in the computer program to represent the values of P along the 
(j — Dth level by the array O[I], I = 1,2,..., NI, and along the jth level by the array 
PID, I = 1,2,..., N1, (where we have assumed N1 mesh points in the x-direction) 
than to use a two-dimensional array. In this notation the implicit scheme becomes 


—Ro* Pisi jer + (1L REY + 47) Pea — Ro” P;-1,j+1 = Rat P[I + 1) 
— R(a* + a )P[I] + Ra PI — 1]+ Ra*O [I + 1] (3) 
+ (1 — Rat + a DOM + Ra O[I — 1], 


where R = 3xr and 
at = 1 + 2e(P[I + 1] + PID), 
a7 = 1 + He(P[T] + PU — 1]. 


We can see from equation (3) that the implicit scheme involves three unknowns 
Pi+1,5+1> Pijer and Pics, and therefore gives rise to a tridiagonal system of 
equations. These can be solved by the recurrence relation form of Gauss elimination 
as given in S: pages 20 to 22. There, the tridiagonal set of equations have the form 


bu — Cit zd 
Ayu, + bu; — C2U3 =d, 
=An-1ly-2 buc. = dy. 


In the computer program we have used the same notation (where the lower case 
letters have naturally had to be changed to upper case). Thus in general 


BI) 2 1 + R(a* +07), 

A[I] = Ra’, 

C(I] = Rat, 
and Df] is given by the right-hand side of equation (3). The solution of the system 
of equations is stored in the array P after the contents of P have been transferred to 


the array O. In this way the solution along the jth level is always in the array P and 
the solution along the (j — 1)th level is always in the array O. 


Since the implicit scheme is a three-level scheme we need values along j = 0 and 
j= 1 before it can be used. Along j = 0 we have initial values already given. For values 


along j 2 1 we use the explicit scheme. However, this process involves some careful 
computation as we shall now illustrate. 


The explicit scheme has a restriction on the size of time step that can be used (in ord 
to retain stability), namely E 


kr <4, rz kh, 
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or 


h? 
T hast 
ks 2k (4) 


The object of using the implicit scheme is that we can use larger time steps than with 
the explicit scheme and thereby reduce the amount of computation required. This 
implies that the time step we use in the explicit scheme will be different (smaller) 
from that used in the implicit scheme. Therefore, to be able to use the explicit scheme 
to "start off" the implicit scheme the step length of the explicit scheme has to be such 
that the step length of the implicit scheme is a multiple of it. The computer program 
therefore allows you to specify any value for the implicit time step (call it kj) and then 
calculates the largest explicit time step (call it kg) such that 


kj=Mkg where M is an integer 
and such that ky satisfies equation (4). 


In fact we choose to use the explicit scheme to provide values along j — k, "and 
j = 2k, (rather than just j = kj). 


This is done to ensure that the effects of the discontinuity in the initial and bound- 
ary data at (1,0), which have been minimized by initially using a finite-difference 
scheme with a small time step, are not allowed to recur at the start of the implicit 
scheme. 
10 Exercises 
The first four exercises discuss aspects of the modelling. 
1 What are the response times for tubes of I in bore which are 

(a) 300 ft long? 

(b) 600 ft long? 

Use the implicit scheme with 7 mesh points in the x-direction and take & = 0.01. 


Hint: to reduce computation time use a step length which requires only about 20 
steps with the implicit scheme. You can estimate the required value by referring 
to the graph of response time against tube length given in the paper by Jones and 
Jordan. 


2 Repeat exercise 1 for tubes of å in bore. 


3 What is the effect on the results of exercises I and 2 if you put £ = 0.1 (the maximum 
value for & in practice)? 


4 Is the linear case (e = 0) a good approximation to the given problem for all 
practical values of £? 


Hint: compare the results of exercises 1, 2 and 3. 
The following exercises discuss aspects of the numerical analysis. 


5 Put L=300, A = 0.1, E = 0, NI = 7, K 2 0.5 and compare the times taken by 
the explicit and implicit schemes to evaluate the response time. 


6 Choose some typical input values for the program and opt for the full print-out. 
Now change statement 670 to 


670 O[NI] - 0 


and run the program with the same data. Compare the ‘ERROR’ columns in the 
two print-outs and say why we choose to take Py,,; = 0.5 when the true initial 
conditions of the problem specify Py, 9 = 0. 


Hint: what is peculiar about the point (I, 0)? 


Note Before going on to other exercises make sure you change statement 670 
back to 


670 O[NI] =.5 
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CONVERGENCE 


This programme studies the circumstances under which the solution of a finite- 
difference scheme cenverges to the analytic solution of the approximated partial 
differential equation as the mesh spacings are reduced. The approaches to convergence 
for initial-boundary and pure boundary value problems are contrasted by investigat- 
ing a simple example from each category. The presenter is Peter Thomas. Reference 
is made to the following notes. 


t y4 

4 

ZA 

Å 

2 AMAA x Z PE 
a initial-boundary value problem pure boundary value problem 
ig. 
Consider the initial-boundary value problem 

ôU OU 0 

å 6d <x<le>0, 

U(x, 0) = f(x) O<x<1, (1) 


U(0,) = U(1,) 20 1>0. 


The partial differential equation can be approximated by the simple explicit finite- 
difference scheme 


Wei m Mg Tui — 24; + Ui+1,) (2) 


where m kjh?, hand k being the mesh spacings in the x- and t-directions respectively. 
The finite-difference approximation to the problem (1) is then : 


Hope m Mgt Muza — 2u;j + Uis) i212..,N— l;j=1,2 
uio = fh) i2L2...,N—1, (3) 
toj = Uyy = O- j=0,12,..., 

where Nh = 1. This can be put in matrix form as 
Uys, = Au; j=01,2..., . 
up =f, (4) 


where u; is the column vector (u; pU Uyiy DÅ 
. " bs "abs TRE 
and Å is a square matrix of order N — I, namely 


1—2r r 


= (SO) fQh),..., SN — DA) 
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The partial differential equation in (1) could alternately be approximated by the 
Crank-Nicolson implicit scheme, 


TM jar b (2 + 200 ay — Uisa gjer 
= ri + (2 — ru + Flisi (5) 
with which the matrix form of the problem becomes 
Buj,; = Cu, J=0,1,2,..., ( 
6) 


uo =f 


where u and fare defined as in the explicit case, B and C are square matrices of order 
N — 1 given by 


[242r =r 7 
=r 2+2r —r 
B= Se 4 Å 
-=r 2+2r —r 
- —r 2+2rl 
[2—2r r 7 
r 2-2r r 
C= $ RS 
r 2—2r r 
L r 2 — 2r | 


If B is non-singular then equations (6) can be solved uniquely and we can write 
uj,; = Du; j20,L2,..., 
au E () 
uo =f 
where D = BC. This is of the same form as equations (4), so that if we drop the 


definition of D as a product of two known matrices B7! and C we can regard (7) as 
representing both the explicit and the implicit method. From (7) we obtain 


Up =f, 

u, = Duo, 

uz = Du, = D(Dug) = D*u,, (8) 
u, = Du,- = D(Du, ;) =... = D"ug. 


Suppose U(x, t) is the true solution of problem (1), and that U; = (U(h, jk), U(2h. jk), 
+++ U[(N — Nå, jk]). Since U(ih, jk) differs from u,; in general, Uj,, — DU; = 
T; # 0, where T, represents the errors incurred by using the finite-difference scheme 
to progress from the jth to the (j + 1)th time level. 


We have 

Ui = DU; +T; j-012.. (9) 
It also follows from (1) that 

U =f. ' (10) 
Writing out (9) explicitly, 

U, = DU, + To 

U, = DU, + T, = D(DU, + To) + T; = D*Uy + DT, + Ty, (11) 


n=l 
U, = D'U, + X DIT, 


j=0 
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Combining the second equation of (7) with (10), 
Up =U =f, (12) 
so that 


aet 
U, —u, = $ DIT, n=1,2,.... (13) 
j=0 
The finite-difference scheme is defined to be convergent to the true solution on the 
fixed time level t if 
lim |[U, — ul = 0, (14) 
h,k =0 
nk-t " 
where || || denotes the norm of the vector (a measure of its size). It follows from (13) 
that convergence is guaranteed if 
nal 
lim} Y D'^!-7T; 


hk 0) j= 
nker || 379 


=0. (15) 


If it takes n steps of length k to reach the fixed time level t (as we have assumed above), 
then 


t = nk. (16) 


Keeping t fixed and letting k ^» 0 means that n ~ oo, so that (15) involves an infinite 
sum. The convergence condition will be satisfied if both ||T;| +0V/ as h,k ~0 
(consistency), and all the powers of the matrix D are bounded (in some sense not 
defined here; it is sufficient for the spectral radius of D to be less than one). À scheme 
for which all powers of D are bounded is said to be stable. Thus consistency and 
stability together imply convergence. This illustrates Lax's Theorem. 


Consider next the pure boundary value problem 


QU OU 
FT ay D<x<1l,0<yx<1I1, 
Ux, 0) = fil)  O«x«l 
U(x, 1) = folx) O<x<1, (17) 
U(0, y) = BO) O<y<l, 
Ul, y) = faly) O<y<l. 
The unit square is covered with a square mesh. The mesh spacing is h = 1/N. 
rt 
Nh=1 
FT 
h 
2 kK —] 
h | 
0 h 2h des Nh=t r 
Fig. 2 
The partial differential equation is replaced by the five-point formula 
Wir + Mang t ijo, + ugs 7 Au; = 0. (18) 
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By writing this down for every internal mesh point we get the matrix equation 
Au = b, (19) 
where u = (tii, Hi 5,... »H3,N—15 12,1, - Uy 1, 1), b has elements determined by 


the boundary values, and A isa square matrix of order (N — 1)? which has the par- 
titioned form 


[B I 
P I B I 
real x 
L I B 
I being the identity matrix of order N — 1 and B the matrix 
[7-4 1 
1-4 1 
B= M . 
1-4 1 
L I -4 


Since all the mesh values are to be computed simultaneously the concept of stability, 
which we used when investigating an initial-boundary value problem, is not meaning- 
ful here. We must therefore approach convergence differently. 


It is shown in Unit I I that when the five-point formula (18) is applied to the Dirichlet 
problem for Laplace’s equation in a rectangle it gives a unique solution. It is also 
shown that if v is any function defined on the set of mesh points in some rectangular 
region D, one of whose sides has length a, then 

max |v] < max |v| + 4a? max |Lv|, (21) 

Då Cy Da 

where D, is the set of mesh points interior to D, C, the mesh points on the boundary 
of D, and L the five-point finite-difference operator (i.e. the left-hand side of (18) 
could be written as Lu; j). 


As the vector u which satisfies equation (19) is only an approximation to the true 
solution U (as in the previous case) we can find a non-zero vector T such that 


AU =b+T. (22) 
Subtracting (19) from this, 
AU — u) = T, (23) 


and since we know the solution of (19) to be unique the matrix A can be inverted, 
giving 
U — u = A`'T. (24) 
For convergence, we require that 
lim |U — ul = 0, (25) 
ino 
and from (24) this will hold if and only if 
lim A^! T] = 0. (26) 
h~o 
Analogously to the initial-boundary value case there are two conditions which when 
satisfied together will ensure that (26) holds. The first is that the elements of T (which 
are in fact the local truncation errors) must tend to zero as h ~ 0. The second is that 


A must be non-singular. The inequality (21) can be used to obtain an upper bound on 
the global error of the finite-difference scheme- This is in fact of the same order of 


magnitude as the local truncation errors. 
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Problem 


Numerical methods can often reveal or illustrate results of mathematical analysis. 
One such result is that a problem requiring the solution of a hyperbolic equation in a 
closed region and the attainment of given values on its boundary is not properly posed. 
Illustrate this by trying to solve numerically the problem 


På a o O<x<1l,0<y<l, 
U(x,0) = fi) O<x<1, 
U(x, 1) = h(x) 0«x«l, 
U(0, y) = fa) O<y<1, 
UL, y) = fio) O<y<l 


Take a square mesh with spacing I and replace the second derivatives with the usual 
central difference expressions. The inverse of the matrix which represents the applica- 
tion of the finite-difference scheme at each interior mesh point must then be found. 
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SHALLOW WATER WAVES 


This programme investigates the behaviour of. standing shallow water waves in a tank 
with vertical ends and constant longitudinal cross-section. The exact period of the 
simplest such wave is found for the case of a tank with horizontal bottom. Results 
from Unit 13 are then used to estimate the corresponding quantity for a tank with 
sloping bottom. 


Introduction 


Ralph Smith opens the programme by displaying water waves in a glass tank. The 
simplest such motions are the stationary or standing waves, which correspond to the 
eigenfunctions of a Sturm-Liouville system derived (later in the programme) from 
the shallow water wave equation. The essential feature of a standing water wave is 
that the surface contains points which neither rise nor fall with time (the nodes). A 
model shows the behaviour of such a surface with just two nodes. 


For the purpose of describing mathematically the water surface motion in a tank with 
vertical ends and constant longitudinal cross-section we take the y-axis to coincide 
with the left-hand end of the tank and the x-axis to lie along the equilibrium level 
of the surface. The configuration of the surface in motion is given by y = n(x, t), and 
the shape of the bottom of the tank (which will not in general be horizontal or even 
flat) is described by y = — h(x). 
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The simplest standing wave has only one node, the surface in this case undergoing a 
slopping motion. We concentrate on this motion (which corresponds to the first 
eigenfunction of the associated Sturm-Liouville system) for the rest of the programme, 
seeking to find its period P. 


The discussion so far has been in terms of the surface shape y(x, t); it is in fact easier 
to deal with the equation for the horizontal fluid velocity u(x,t) and this is the 


dependent variable used hereafter. u and y are connected by the two equations 
U= —8Nx, 
(1) 
(hue = =n, 


where g is the acceleration due to gravity and subscripts are once more being used to 
denote partial differentiation with respect to the subscripted variable (see Appendix 
for a derivation of these equations). 


The horizontal fluid motion at the surface of the water in the studio tank is demon- 
strated with the aid of markers. 


The Shallow Water Wave Equation 


Dominic Jordan introduces the shallow water wave equation 
1 
(hu)... = ge (2) 


which is obtained by eliminating y between the two equations (1) above. This equation 
describes all possible motions within the limitations of the model. Since h is not in 
general a constant the equation is more complicated than the wave equation for (say) 
a uniform string. The horizontal fluid velocity is clearly zero at the ends of the tank, 
so if | is the tank’s length we have the boundary conditions 


u(0, t) = ull, t) = 0. (3) 


Ioue]. gud), 
uo, =ull,)=0 
Equations (2) and (3) together lead to an eigenvalue problem if we try and solve by 
separation of variables. Putting 
u(x, t) = X(x)T(t) (4) 


effectively restricts our attention to the standing wave solutions, as is demonstrated 
by a piece of film. Solutions of this form always have their nodes, maxima and minima 
(with respect to variation of x) at fixed places, keeping the same basic shape at all 
times; only the amplitude varies. 
Substitution of (4) into (2) and (3) gives 

(hX)" + AX = 0, 
(5 
X(0) = X(I) = 0, 
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together with . 
T" + ÀgT — 0, (6) 


where å is the separation constant. 


Tank With Horizontal Bottom 


We consider first the case where h(x) = H, a constant. This corresponds to a tank with 
horizontal bottom. Equations (5) can then be written as 


À 
bap ege 
X+ 0, 


X(0) = X(I) = 0. 


Ralph Smith points out that the solutions (i.e. eigenfunctions) of this system corre- 
spond to standing waves in the tank. The eigenvalues of the system are 


(1) 


2,=ns?H/P  n=1,2,..., (8) 
and the corresponding eigenfunctions are 

X,(x) = sinnmxl)  n=1,2,... (9) 
Solving equation (6) then gives 

T(t) = C, sin [(nnH*g*t/I) + c,] Hn-1,2,..., (10) 
so that 

u(x, t) = C, sin [(nzH?g?t/I) + c,] sin (nax/)  n=1,2,..., (11) 


is a set of separated solutions of equation (2) satisfying the boundary conditions (3). 
C, and c, are arbitrary constants. C, is the amplitude of the motion at a maximum 
or minimum of X,, and c, depends upon what stage of the cycle has been reached at 
tsu 


Xiah X=0 

XIX) 0 

T*4T-0 a 
procu 3 
ix re Cainfnx Hg iD cnl 


These constants would normally be determined by initial conditions, but typical 
initial conditions will not produce standing wave solutions. Since we are interested 
mainly in the periods of the motions represented by (11) there is no need to specify 
initial conditions. 


We get the simplest case from (11) by putting n — 1. This gives (dropping the subscript 
throughout) s 


u(x, t) = C sin [(xH*g*t/I) + c] sin (zx/D. (12) 


A film shows the motion of the velocity profile with time in this case and the corre- 
sponding behaviour of the water surface (the simple slopping motion described 
carlier). It follows from (1) that the shape of the surface is given by 


n(x, t) = C(H/g)* cos [(nH*g*t/I) + c] cos (nx/I) (13) 
(see Appendix for details). The period of the motion is seen to be 
2 21 
2 fH (14) 


nHigP Hig? 


Assigning typical values I = 1m, H = 0.04m, g ~ 10ms~?, we find P to be somewhere 
in the region of three seconds. 
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Tank With Sloping Bottom 


Dominic Jordan considers the more complex situation which arises if the bottom of 
the tank is not horizontal. The equilibrium depth is taken to vary linearly from 2cm 
at one end to 4cm at the other, and it is supposed that the tank is Im long. The depth 
function is therefore 


h(x) = 0.02(1 + x) O<x<1. (15) 


Fe 


A(x)=0.02(1+x) ,O<x<1 


The situation is modelled by equations (5) (with | = 1) and (6) together with (15), 
assuming once again that attention is being restricted to the standing wave solutions. 
In order to tackle (5) using the methods of Unit 13 the differential equation must be 
put into self-adjoint form, i.e. into the form 


(pX'Y — qX + ApX = 0. (16) 
We find p = h?, q = — hh", p = h, so that the differential equation becomes 
(XY + hh X + AhX — 0. (17) 


Substitution for A from (15) gives the Sturm-Liouville system 
(Ll -xPX]--A5001-xX-0 O<x<1, 


18 
X(0) = X(1) = 0. us 


[Note that in the notation of (16), p = (1 + x)?, q = 0 and p = 50(1 + x), so that the 
form of the differential equation in (18) satisfies the conditions p>0,q20,p>0 
assumed on W160. We can therefore be sure that the theorems of Unit 13 will apply 
to this system.] 1 


The problem (18) cannot be solved exactly, but something can still be said about the 
shape of its eigenfunctions. 


po 


COS DENDUM 
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We know from the Oscillation Theorem [rather than from the Separation Theorem 
as is stated] that the kth eigenfunction of the system (18) has exactly k — 1 zeros in 
(0, 1). This enables us to sketch the first four eigenfunctions as shown in the left-hand 
half of the photograph above. [Clearly the sketch for X,, is in error: we require X ,(1) = 
0.] The pattern continues for later eigenfunctions, Using equations (1) we can deduce 
the corresponding water surface configurations n,(x, to), n = 1,2,..., at some fixed 
time t = ty for which (x, t) is not momentarily passing through its equilibrium 
position. The zeros of X, correspond to the maxima or minima of TX; to). [This 
follows from the first of equations (1). We see from the second equation that the zeros 
of 71,(x, to) correspond to the maxima or minima of (1 + x)X,. It is not therefore true 
as stated that the zeros of r,(x, to) correspond to the turning points of X,. There must 
however be a turning point of (1 + x)X, between any pair of consecutive zeros of X,, 
so we may justifiably conclude that n,(x, ta) has a zero between any two consecutive 
zeros of X,,, See Appendix for further details.] We obtain the shapes of water surfaces 
corresponding to the first four eigenfunctions of (18) as shown in the right-hand half 
of the photograph; the slope of 7,(x, tọ) is zero at cach end of the tank. 


Knowledge of the eigenvalues allows us to calculate the period of the corresponding 
standing wave motions. Equation (6) has the solution 


T(t) = C sin [(Ag)*t + c], (19) 
and this has period 
P = nf. / Ag. (20) 


The period of the standing wave corresponding to the nth eigenfunction of system 
(18) is therefore found by substituting the nth eigenvalue for 7 in this equation. We 
now concentrate on estimating the lowest eigenvalue 2, which will provide a corre- 
sponding approximation to the period of the simplest standing wave. 


The Monotonicity Theorem is used to place upper and lower bounds on År, We see 
that in the interval [0, 1] the functions p(x) = (1 + x)? and p(x) = 50(1 + x) satisfy 
the inequalities 
1 < p(x) <4, 
P(x) 21) 
100 > p(x) > 50. 


The theorem then tells us that the nth eigenvalue of the system (18) is greater than the 

nth eigenvalue of (16) with p = 1, q = 0, p = 100, and less than the nth eigenvalue of 

(16) with p 2 4,q = 0, p = 50, the same boundary conditions pertaining in each case. 

The first eigenvalue of (18) therefore lies between the first eigenvalues of the systems 
X” + 1004X 0  O<x<1, 


Q2) 
X(0) = X(1) = 0, 
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and 


4X” + 504X = 0 O<x<1, 


(23) 

X(0) = X(1) = 0. 
These first eigenvalues are readily found to be n?/100 and 472/50 respectively, so that 

= LAS z (24) 
or 

0.098 < 1, < 0.790. (25) 
Substituting these bounds into equation (20) we find 

6.39s > P > 2.255. (26) 


This is hardly a spectacularly accurate estimate for the period but it’s only a first 
attempt. 


Conclusion 


Ralph Smith states that the procedure needed to improve the bounds found for the 
period P of the simplest standing wave is one of trial and error. He previews Radio 
Programme 7 which is devoted to finding better bounds on P. The problem is 
approached using the minimum principle 


[ wor + att 

A, ——;, (27) 

[oos 
0 


where p(x) = (1 + x}. q(x) = 0, p(x) = 50(1 + x), and $ is any continuous piecewise 
continuously differentiable function which vanishes at x = 0 and x = 1. The par- 
ticular function ¢ used in the radio programme is 


$é(x)-xl—-x O<x<l. (28) 


APPENDIX 
Derivation of the Shallow Water Wave Equation 


Consider a tank with constant longitudinal cross-section containing liquid. As before 
we take the x-axis to lie along the equilibrium level of the liquid. The surface profile 
of the liquid is given by y = (x, t) and the shape of the bottom of the tank by 
y = h(x). 


We make the following initial assumptions: 


(a) The motion of any particular fluid particle lies at all times in a plane parallel to 
that formed by the x- and y-axes, i.e. we consider the fluid motion to be two- 
dimensional. 
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(b) The fluid is incompressible, so that its density p is constant. 
(c) The effects of viscosity are negligible. 


(d) The vertical component of acceleration of the liquid particles has a negligible 
effect on the pressure; in other words the expression for the pressure p(x, y, t) at 
the point (x, y) of the fluid at time t is taken to depend linearly on depth as in the 
hydrostatic case, so that 


p(x, y, t) = pgln(x. 1) — y], (29) 
where g is the acceleration due to gravity. 


(e) There exists at least one value of t for which the horizontal component u of the 
liquid partiele velocity does not vary with depth. 


Assumption (d) is basic to shallow water theory. Differentiating equation (29) with 
respect to x, 


Px = pens ^ (30) 


from which we see that py, the pressure gradient along the tank, is independent of y. 
This pressure gradient is however (up to sign) the x-component of the force per unit 
volume, so we may deduce from Newton's Second Law that the horizontal component 
of particle acceleration, du/dt, is independent of depth. It follows from assumption (e) 
that u is independent of y at all times. This means that particles lying in any plane 
perpendicular to the x-axis will remain in a moving vertical plane as they sweep to 
and fro. 


Consider the fluid flowing into the space between two fixed vertical planes P through 
x and P' through x + Ax. We can form an equation of conservation of mass by 
equating the net inflow rate of fluid into the region between the two planes with the 
rate of volume increase due to the vertical velocity of the surface. 


The volume flow rate from left to right through P is b[n(x, t) + h(x)]u(x, t), where b 
is the width of the tank. The volume flow rate from left to right through P' is 
b[n(x + Ax, t) + h(x + Ax)Ju(x + Ax, t), so that if Ax is small the net rate of inflow 
into the region between P and P' will be ` 


b[n(x, t) + h(x)Ju(x, t) — b(n(x + Ax, t) + h(x + Ax)]u(x + Ax, t) 
œ —b[(n + hjul Ax. (31) 


The vertical velocity of the surface at x is 7,, and the corresponding rate of change of 
volume between P and P' is approximately 7,bAx. Equating this with the right-hand 
side of (31), 


[o + hu), = =r (32) 


We obtain a second equation relating u and y by writing down Newton's Second Law 
for motion in the x-direction. The particle acceleration in this direction is du/dt = 
u, + uu, (compare with line I on page 15 of Unit /) and the force per unit volume is 
— px. Hence 


p(u, + un) = — Py- (33) 
Using equation (30) this becomes 

u, cru. = —gny. (34) 
We have two equations connecting u and y. The additional assumption is now made 
that u and y. together with their derivatives, are small quantities whose squares and 
products can be ignored in comparison with linear terms. When this simplification 
is made equations (32) and (34) yield equations (1) of these Broadcast Notes, 

U= gl (1) 

(hu), = —n- 
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Differentiating the first of these with respect to x, the second with respect to t, and 
equating the two expressions for —1,, gives the shallow water wave equation 


1 
(hu), = -u,. Q) 
8 


A more comprehensive derivation of equations (1) can be found on pages 23, 24 of 
Water Waves by J. J. Stoker (Interscience, 1957). 


Derivation of n,(x, t) From u,(x, t) = X,(x) T. (t) 


It seems intuitively reasonable that the function ,(x, t) corresponding to the standing 
wave solution u,(x, t) = X,(x)T,(t) should itself represent a standing wave (and hence 
be a product of a function of x and a function of 2). It is not however mathematically 
obvious that this need be so, and it is clear that such a function would not be the 
unique solution to equations (1). 


If 7) and n(” are assumed to be two distinct solutions of (1) then 7, = nP — 1 
must satisfy (Fa) = (fi), = 0, whence any two solutions may differ by a constant. We 
therefore proceed to construct a separated solution for rj, (x, t) on the assumption that 
there is one, and then show that the constant which may be added to this solution to 
form a further solution of (1) must in fact be zero from other considerations. 


Suppose u(x, t) = X,(x)T,(t), where X, is the nth eigenfunction of the system (5) and 
T, the corresponding solution of equation (6). Suppose further that 7,(x, t) = F,(x)G,(0), 
and that equations (1) are satisfied by u = u,,4 = 1. We have 


(X, T, = —F,G,, 


35) 
X = -gEG, 
Separation constants o,, f, can therefore be found such that 
QGXy Ga 
Roc. ^ 
3 
E R (36) 
m uw 
and these in turn give 
F, = (hX p) lans 
7 (37) 
G, = Tifgf,. 
so that 
n, = XY Tr fgasf),. G8) 


Differentiating the first of equations (37) and substituting for F; in the second of 
equations (36) we see that 


(hX) + of, X, = 0, (39) 
and it follows from (5) that PN = i,. Hence 
tn = (AX) Ta Ang (40) 


As has been mentioned a constant may be added to the right-hand side of (40) and the 
resulting expression will still satisfy equations (1). Suppose ¥,(x, t) = nx, t) + c. The 
volume of liquid in the tank must always be the same and consequently equal to the 
volume in the equilibrium state. For all t therefore, 


l 
Í nx, Ddx = 0 (41) 


if}, is to be a physically satisfactory expression for the surface profile. Since 


I 
Í ital, t)dx = [AX — h(0)X,(0)]7,(0/2,8, (42) 
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Which vanishes by the boundary conditions of the System (5), we see that c = 0 and 
that 5, = 7. Equation (40) therefore gives the unique expression for n,(x, t) corre- 
sponding to u(x, t) = X,(x)T;(t). In particular it gives (13) as the surface profile 
corresponding to the velocity function (12), since in this case 4, = z?H/I?. 


To Show That the Zeros of u, (x, to) and 1,(x, to) Interleave 


We wish to show that exactly one zero of nx, to) lies between each pair of consecutive 
zeros of u,(x, to). From equation (40) this is equivalent to showing that exactly one 
zero of (hX „y lies between each pair of adjacent zeros of X,,. Since all the zeros of X, 
are zeros of hX, there must be at least one turning value of hX, between any two 
consecutive zeros of X,,, i.e. there must be at least one zero of (hX,y. 


X, must be of constant sign between any two of its adjacent zeros. Suppose without 
loss of generality that X, > 0 between two such zeros. It follows from (5) together 
with the inequality 2, > 0 that 


(hX,)" <0 (43) 


between the two zeros in question. Any turning point of hX, in this interval must 
therefore be å maximum, and since two maxima of a continuous function are neces- 
sarily separated by a minimum there can only be one such maximum. There can 
therefore only be one zero of (hX,) between the two zeros of X,. 


Note incidentally that the inequality (43) and the reverse inequality for the case when 
X, < 0 provide more information with which to sketch the graphs of the X,,. 
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RADIO PROGRAMME 7 Broadcast Notes 


SHALLOW WATER WAVES 


The investigation of Television Programme 3 into shallow water waves is continued. 
After a recapitulation of material in TV3, improvements on the bounds previously 
calculated for the period of the simplest standing wave motion in a rectangular tank 
with sloping bottom are sought. Successively more restrictive bounds are found by 


(a) applying the minimum principle to the first eigenvalue of the Sturm-Liouville 
system which arises from the shallow water wave equation, 


: (b) looking at a different Sturm-Liouville system which has the same eigenvalues as 
the system in (a), and bounding the first eigenvalue using the Monotonicity 
Theorem. 


(c) applying the minimum principle to the first eigenvalue of the system in (b). 


In each case the bounds on the period can readily be calculated from those on the 
first eigenvalue. 


The presenters are Ralph Smith and Dominic Jordan, and the introduction is by 
Ralph Smith. Reference is made to the following notes in the programme. 


h(x)=0.02(1+x), O<x<1 


Fig. I 


Putting u(x, t) = X(x)T(t) into the shallow water wave equation 
" 1 
(tek = Up, 
E 


where A(x) is given in this case by the expression in the photograph above, and using 
the boundary conditions u(0, t) = u(1, t) = 0, we find that 


[(1 + x)X]" + 150X = 0, (1) 
X(0) = X(1) = 0, (2) 
T"  ÀgT — Q. (3) 


After equation (1) has been put into self-adjoint form we have 


[A + xX] + 25001 + x)X = 0, (4) 
X(0) = X(1) = 0, (5) 
T" + ÀgT — 0, (6) 
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and the first two of these equations form a Sturm-Liouville system of the type con- 
sidered in Chapter VII of W. Corresponding to each eigenvalue 2, there is an eigen- 
function X,(x) and a standing wave u,(X, t) = X,(x)T,(t). The period of this standing 
wave is the period of 7,(), which is a solution of equation (6) found when 4 = 4, 
namely 


T) = C sin [G/2,g) + c). am 


[The constants C and c should have subscript n in this context but it has been 
omitted.] This has period 


Py = In f Ag. (8) 
The complete standing wave corresponding to 4, has the form 
ux, 1) = Xx). C sin [G/A,gt) + c] 9) 


We now concentrate on the first eigenvalue 4,. This corresponds to the simplest 
Standing wave u(x, t) which has period P, = 2n] f Ag. It was found in Television 
Programme 3 that 


n? 


0.098 < A, < 0.790, (10) 
and consequently that 
6.39s > P, > 2.25s (11) 


(see inequalities (25), (26) of TV3 Broadcast Notes). 


We next try and improve these bounds using the minimum principle 
1 
NES 
hg =o 
| gax 


o 


(12) 


where in this case p(x) = (1 + x)?, q(x) = 0 and p(x) = 50(1 + x). $ is any continuous 
piecewise continuously differentiable function which vanishes at x = 0 and x = 1. 
Equality holds if and only if ¢ = X,. We might suppose that (x) = sin zx would 
give a close approximation to Å, since this was found to be the first eigenfunction in 
the case of constant depth (put I = n = 1in equation (9) of the TV3 Broadcast Notes). 
The integration involved is a little lengthy however so we choose instead the function 


Ox) =x(L-x)  O<x<1, (13) 


which has a similar graph to that of sin nx but ensures that all the integrands in (12) 
will be polynomials. 


We get Fig. 2 
1 
| (1 + xy1 — 2x)*dx 
4, <p, (14) 
f 50(1 + x)(x — x?)?dx 
0 
which gives 
À, < 0.320. (15) 
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Putting this together with (10) we see that 
0.098 < 4, < 0.320, 
and hence 
6.39s > P, > 3.558. (16) 


[Putting #(x) = sin x produces an upper bound of (147? + 3)/450 = 0.314 for Å 
and a corresponding lower bound of 3.58s for P,, so hardly any improvement is 
gained using this expression for Ø(x) instead of that in equation (13).) 


The problem can be reformulated by putting 


Y — (1 3)X, ; (17) 
and substituting into equations (1) and (2): 
504 
r" Y=0 2 18 
y+ Tox ^ (18) 
Y(0) = Y(1) = 0. (19) 


The eigenvalues of the system (18), (19) are the same as those of the system (1), (2), for 
if X,(x) is the nth eigenfunction of (1), (2) it will satisfy equation (1) with A = 4,. 
Y(x) = (1 + x)X (x) will then satisfy (18) with A = À,,and Y(0) = X,(0) = 0, ¥,(1) = 
2X,(1) = 0. Equation (18) is moreover already in self-adjoint form. Equation (3) 
remains unchanged by the transformation (17): 


T" + AgT =0. (20) 


Comparing equation (18) with the standard form 
(pX'Y — qX + ApX = 0, (21) 


we see that in this case p(x) = 1, q(x) = Oand p(x) = 50/(1 + x). Since 25 < p(x) < 50 
on the interval [0, 1], we know by the Monotonicity Theorem that 4, must lie between 
the first eigenvalues of tHe systems formed by adding to the boundary conditions (19) 
the differential equations 


Y" + 50AY = 0, (22) 

Y" + 252Y = 0. (23) 
The first eigenvalues of these systems are 27/50 and n?/25 respectively, so that 

0.197 < A, < 0.395. (24) 


Combining this with our previous best estimate (the inequalities prior to (16)) we have 
0.197 < 24 < 0.320, (25) 
or 
4.528 > P, > 3.55s. 


The minimum principle for the Sturm-Liouville system (18), (19) is 


ÅS 35% (26) 


and using f(x) = x(1 — x) once more, 


1 
f (1 — 2x)*dx 


0 
hi SKE (27) 
Í —————dx 
o 1+x 
which works out to give 
Ay < 0.295. (28) 
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Comparing this with (25), 
0.197 < 4, < 0.295, 
so that 


4.525 > P, > 3.695. 


Another way of estimating 4, is suggested by SAQ 23(c) of Unit 13, where it is shown 
that 


1 2 


lim a = 7, (29) 
ES | f der] 


where 4, is the kth eigenfunction of the problem 
u” + Apu = 0 x «x «f, 
u(a) = u(B) = 0. 


The expression on the right-hand side of (29) is a reasonable estimate for Auk? even 
when k is not large provided that p is not too irregular. Putting æ = 0,8 = 1 and 
(from equation (18)) p(x) = 50(1 + x)^!, we arrive at the value 0.288 for the right- 
hand side of (29), which may be taken as an approximation to 4, 


It is in fact possible to find 4, explicitly using the methods of Unit 14. It takes the 
approximate value 0.290, so that the upper bound of (28) and the value derived from 
(29) are seen to be quite good estimates. The corresponding value of P, is 3.73 s. This 
compares with values for a flat-bottomed tank of length 1 m and water depths 2 cm, 
3 cm, 4 cm of 4.57 s, 3.69 s and 3.20 s respectively. 
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TELEVISION PROGRAMME 4 


REVISION 


This programme could loosely be described as adopting an algorithmic approach to 
the solution ofa partial differential equation together with its associated initial and/or 
boundary conditions. The intention is to review the different types of equation which 
have arisen during the course, to look at the allowable initial/boundary conditions 
associated with each type, and to enumerate the possible analytic and numerical 
methods of solution available in each case. The programme is introduced by Daniel 
Lunn, and the second presenter is Peter Thomas. 


Categorizations 


The first part of the programme makes periodic reference to the equation 
Uys — Uy = SEU (1) 


(where subscripts again denote partial differentiation with respect to the subscripted 
variable), and after a solution domain and certain initial conditions have also been 
specified the resulting problem is solved. 


An initial and/or boundary value problem consists of both a partial differential 
equation, whose solution is required in some given domain, and a set of conditions 
to be satisfied by the solution on the boundary of this domain. Three broad categories 
of initial/boundary conditions are considered in the course, namely pure initial 
conditions, mixed initial-boundary conditions and pure boundary conditions ; these 
may be symbolically represented as shown in the photograph. [The situation at top 
right is stated to represent a fourth distinct type of condition and a second kind of 
pure initial condition. Initial conditions of this type have not really been investigated 
in the course and do not in any case difler substantially from those represented by 
the diagram at top left.] The boundaries in the cases represented by the two lower 
diagrams might be at infinity, in which case the associated conditions would be 
expected to hold in the appropriate limits. 
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As a first step towards solving equation (1) we seek to classify it as hyperbolic, 
parabolic or elliptic, and then look for a method of solution. The three principal 
analytic methods mentioned in the course are those of characteristics, Separation of 
variables, and Green's functions (other non-numerical techniques include the finite 
Fourier transform but that is not considered here). 


We could also attempt a numerical solution if an analytic one is not forthcoming. 
Any information gleaned from the analysis, such as the likely effect ofany discontinuity 
in the initial or boundary data, will be helpful when we try to solve the problem 
numerically. We may choose to use either the numerical method of characteristics 
or a finite-difference scheme, but having made this choice we must still decide how 
best to use them. 


A solution domain and initial conditions are added to equation (1) in order to form 
a properly posed problem: 


Uyy — Un = U, + Uy xeR, t>0, 
u(x, 0) = xe#* xeR, (2) 
u(x, 0) = (1 — 4x) xeR. 


Equation (1) is classified by investigating the sign of the quantity B? — 4AC, where 
A, B, C are the coefficients of the Uy, Vy and u,, terms respectively. Since in this case 
A= —1, B = 0 and C = 1, we have B? — 4AC = 4 > 0, so that the equation is 
hyperbolic. 


Hyperbolic 


We know that in general hyperbolic equations can occur in properly posed problems 
of initial condition or initial-boundary condition type, but not in pure boundary 
value problems. We can choose between two possible analytic methods of solution, 
the method of characteristics or separation of variables. The initial and/or boundary 
conditions usually determine which choice to make in order to derive the analytic 
solution in its simplest form. Characteristics are most suitable in pure initial value 
problems (we could also use the Fourier transform method introduced in M201 31, 
but that has not been considered in this course). In initial-boundary value problems 
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both characteristics and separation of variables are applicable but the latter method 
is usually preferred since it tends to reduce the problem to a Sturm-Liouville system. 
A solution to the problem may then be constructed by using the eigenfunctions of 
such a system. 


The problem given by equations (2) is a pure initial value problem, so the method of 
characteristics is used to solve it. The viewer is expected to supply the details of the 
subsequent calculation, but a rough outline is given. The characteristic coordinates 
in this case are found to be 


f=x+tn=x—6 (3) 
and after transformation into the C, 4 coordinates equation (1) becomes 


May = s. (4) 


Pe 


[Io HH. 


The general solution of this equation is 
u(5 n) = f (e? + gn) 6) 


where f, g are arbitrary differentiable functions. [The general solution of equation (4) 
isgiveninthe programme by theexpression shown in the photograph above. Seasoned 
M321 students will have no difficulty in identifying the deliberate mistakes!] 


Transferring back to x, t coordinates, 


u(x, t) = f(x + Det + g(x — t). (6) 
In order to satisfy the initial conditions given in (2) we must choose f(¢) = č + c, 
gin) = — ce, where c is an arbitrary real number. Hence the unique solution to the 
problem is 

u(x, t) = (x + thet, (7) 


The problem (2) would not have been so simple to solve if an additional term, +u, 
had been added to the right-hand side of the partial differential equation. A term 
+4u would appear on the right-hand side of the canonical form (4) [the factor 4 is 
omitted in the programme], and it is not clear how the resulting equation can be 
solved analytically. 


It can however be solved numerically, using either the method of characteristics or a 
finite-difference scheme. The first of these appears to be the better option in this case 
since the equations of the characteristics are already known and are simple. In 
general this method is to be preferred because it closely reflects a known analytic 
property of hyperbolic equations, namely that the solution is propagated along the 
characteristics or is otherwise governed by their geometrical disposition ; in particular 
discontinuities in the initial or boundary data are propagated along the characteristics. 
Discontinuities in the solution cannot be so easily accommodated when using a finite- 
difference scheme. The numerical method of characteristics is described in section 2.3 
of Unit 2. 
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Parabolic 


We next consider parabolic equations, typified by the diffusion equation u, = uyy. 
These may occur in pure initial value problems or initial-boundary value problems 
[although the programme only acknowledges the second possibility] and the usual 
analytic method of solution is separation of variables [the finite Fourier transform 
may be applied in some cases, and this is essentially a Green's function technique; 
the method of M20131 would again be used for a pure initial value problem]. 
Separation of variables leads to two ordinary differential equations, one of which 
forms part of a system whose eigenvalues are to be found exactly or else estimated 
using the methods of Unit /3. 


The only available numerical method for initial-boundary value problems involving 
parabolic equations is that of finite differences. The partial differential equation and 
associated conditions are replaced by finite-diflerence equations and the solution 
‘domain is covered with a rectangular mesh. The solution of the finite-difference scheme 
at the intersections of the mesh is sought. It is calculated on cach successive time level 
using values previously computed for the preceding level, the process being started 
off with a set of known values on the initial line t — 0. Since the finite-difference 
scheme is only an approximation to the original problem, small errors (local trunca- 
tion errors) are introduced into the calculation at each mesh point. We wish our 
Scheme to be stable so that these errors will not accumulate. We also require that the 
accuracy of the finite-difference solution should improve as the mesh spacings hand k 
are decreased, and that the solution should tend to the true solution of the partial 
differential equation as h,k -« 0 (which is the definition of convergence). For this 
to happen the local truncation errors must tend to zero as h,k 0, i.e. the scheme 
must be consistent. In fact we know that, for the class of equations considered in the 
course, stability and consistency together imply convergence, so only the first two 
of these properties need be established. 


Elliptic 


There remain problems involving elliptic equations, such as Poisson's or Laplace's 
equations. In order to be properly posed such problems must be of pure boundary 
condition type. The possible methods of solution include separation of variables and 
Green's function techniques. The former is used wherever possible because it is 
manipulatively easier. Laplace's equation u,, + u,, = 0, being homogeneous, 
separates straight away, but Poisson's equation Uxx + Uy, = p usually presents a 
tougher problem, since the p on the right-hand side could be any function of x and y. 
Even if it is not possible to separate Poisson’s equation directly it may be that 
separation can take place after transformation to new coordinates, or after changing 
the dependent variable. For example the equation 


Ux, + Uy = 1 (8) 
is not immediately separable, but if we define v(x, y) by the equation 

u(x, y) = v(x, y) + På, | 0) 
we find that 

Uy + Vy = 0, (10) 


which is Laplace’s equation once more. 


If no transformation will enable us to use separation of variables we can fall back on 
Green’s function techniques. If Poisson’s equation 


Ux. + Uyy = p (x, y)e Å (11) 


is part of a pure boundary value problem, the Green's function for the system, 
G(x, y;6,n), will satisfy the equation 


G,, + Gy =0 6 y) 7- (Gm (12) 
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` together with appropriate conditions on the boundary of A. The solution to the 
problem will be given as 


u(x, y) = -ff G(x, y :č, n)p(č, n)dč dn. (13) 
A 


[This equation appears in the programme without a minus sign. This is a matter of 
convention, but the convention adopted in M321 (i.e. the precise way in which G 
has been defined there) necessitates the presence of a minus sign.] It is often difficult 
to find G, and even if G is found the integral in (13) may defy attempts to work it out. 
In such cases we might turn once more to numerical methods. 


The finite-difference method is again employed for problems containing elliptic 
equations. The region covered by the mesh is closed on this occasion, and the partial 
differential equation and associated conditions are therefore replaced by the algebraic 
problem of solving a set of simultaneous equations, one for each mesh point. We have 
to consider whether the solution of the algebraic problem approximates the solution 
of the partial differential equation ; we also need to know how to solve the algebraic 
problem. The first of these points is satisfactorily resolved if the solution of the 
algebraic problem tends to the solution of the original boundary value problem as the 
mesh spacings tend to zero (which is convergence once again). The only method for 
Solving the algebraic problem which has been discussed in the course is simple 
iteration (the Gauss-Seidel or SOR methods for example). 
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COURSE NOTES Course Notes 


ERRATA IN SET BOOKS 


The following errata have been observed in the set books. These lists include reference 
to sections not contained in course reading passages. 


Errata in W 


W2, lines — 5, — 3: the expression for the force should read 


TG, 062.0) Te, 926.0 


[is] +(e) 1. [V8] * i] ae 


W4, line 13: delete the full stop 
W14, line — 14 : should read 


u(0, ) = 3L f (ct) - f(- c0] + all g(x) dx = 0 


W26, line I : should read 


aff. jt Sext dé dij 


W27, line 8 : should read 


2r rå 
«sf li (—tsin?zx) dx +--- 
2Jal -3+ 


W33, line 2: for v!, v? read vo, v? 7 
W34. line —15: for L, [f= f, read L,[8] = J, 
W35, exercise 5: should read 


Show that the problem 
Ou ,0u 


3270 for0<x<1lt>0, 
u(0, t) = fi(t) 
Ån = fit) 


u(x,0) = f(x) forü&x&! 
—(x,0) = f(x) forO<x< 


has at most one solution. 


W36, line —5: for fy, fo, read f. g, 
W53, line —13: The third integral in this equation should read 


+ f |ørad ul? dx dy 
D 
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W53, line — 11: for membrance read membrane 

W55, line 9: for somehwere read somcwhere 

W74, line 16: for f = 0 is the only function read functions f which are zero at all but 
a finite number of points are the only functions 

W74, line 17 : for its read their 

W104, line —2: for U read u 

W126, exercise 4: for f(x) read — f(x) 

W133, line 3: delete uniformly in $ 

W138, line 1: for 0 < x < Aread0 < y <A 


m pn A pr 
W138, line 8: for Í f read Í f 
o Jo o Jo 


W142, line 15: should read 


f f Fx, y dx dy = zf agdy + 2] [a2 + bl] dy 
mg e 1 += 
W167, line —6: for right-hand read left-hand 
W172, line —6: fora < x$ åreada & x <% 
W174, line 1: for qb?dx) read qb?) dx 
W175, line 9: It is stated that "We can prove a monotonicity theorem like that given 
above for this problem”. This statement is misleading. Monotonicity holds with respect 
to æ p, q, p and b, but not with respect to Å. 
W179, line 16: insert ) before = 
W180, line 4: for (41.1) read (40.1) 
W182, line 11: for sin nO read sin m0 
W186, equation (42.2): for cy, (t) read c," (t) 
W186, line — 1: for Cno read Cpo 
W190, line — 1: should read 

(=2)" ca kita 1 k—A 

CP ^ (4*1 ^ k+1-4 k*i-À 

W191, line 2: for nc, read n c, 
W194, equation (44.3): for p read à and for sin 0 d0 dọ read sin 0 dd d$ 
W377, line —3: insert ) after k 
W385, line 3: for (83.2) read (82.3) 
W387, line 20: for [sin X + sinX Q read [sin X + sin?x( 


W399, line 5: for oe d read 3x)? 
Ox s 


W402, solution to section 7, exercise 2: should read 


2. The part of the strip 0 < x<}, where 
0<t €2-— llog cos x — log cos 4. 


W402, solution to section 7, exercise 3, line 2: should read 
0<: €3— |an^!x —tan^! 4, 

W402, solution to section 7, exercise 4, line 2: should read 
t > |tan^!x — tan” HI. 


W404, solution to section 10, exercise 2, line 2: for /a/Dz read ./a/Dz 


W404, solution to section 10, exercise 2, lines 2 to 5: for 
If the transformation £;...and AD > 0 
read 
If the transformation č; — ex + fly + y;iz, i= 1,2,3, reduces the three- 
dimensional operator to standard form, choose unit vectors j and v which are 
perpendicular to each other and to y, and show that the transformation £s 


-ax + p- By, Ča =v-ax + v.fly reduces the two-dimensional operator to 
standard form. 


W404, solution to section 10, exercise 3: for V2u read aV?u 
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Wa 17, solution to section 30, exercise 1, line 3: for log (r? read log [r? 
W418, solution to section 32, exercise 3, line 2: should read 


sin(n — $)x sinQm — Dy sinh/(( — 3)? + Qm — 1} + 1] (1 — z) 
Qm — DVI(n — 3)? + Qm— 1? + 1] cosh/[(n — $) + Qm — 1? + 1] 
sin cJ/ Ave sin ev/ At 


YTE 


W423, line — 12: for a sinh č, sin n read o sinh č sin 
WA423, lines —7 to —1: should read 


W422, solution to section 41, exercise 2, line 2: for 


Y” + [A — 00054] Y=0 for0 <y < 27, 
Y(2z) — Y(0)=0 
Y'(27) — Y'(0) = 0 


X"  [x?e?cosh?t — AX 20. for0<é< ann (Å, 
a 
à -i[b å 
X'(0) = X|tanh p =0 ifYQz — n) = Y(n), 


b , 
X(0) 2 de) =0 if Yn — n) = — Y(n). 
The 7-problem can be solved to give eigenvalues 4, (c), A2(w), . . . as functions of c». 


Then the ¢-problem is used to determine a discrete set of values of œ. 
W424, solution to section 44, exercise 4 : should read 


$c 
4. u(r,0, 0,1) = m 2 Pr a Jaami [yer Oy 
n=0 K= 


E 
[oo eos O+ Y (acos mó + b,sin mo)P? (cos o} 
m=1 


1/2 +(1/2 1 +(1/2: ET 
where yup * 0/729 Jim up er ») = Dns DV ue Gry) = 0, 


k= 123 
Qn px pt 
(2n + 1)(n — m)! Í I Í S(r, 0, d)r Jara VE yy Pa(cos 0)r?sin 0 cos mødr dO dø 
o Jo Jo 


Ankm I 
In(n + m [ Ju Vy dr 
0 


W440, lines — 2, — 1: should read 
lu(t, 1) — url, 1)] € 4(cosh 2 — 3) 


= 0.381 
Errata in S 
: ax dx 
S10, lines 4, 5: for ax read ax 


S10, line —7: for (1.7) and (1.8) read (1.8) and (1.10) 
S18, line —9: for (0 <x € I) read (0 < x « I) 
S26, line —4: for ut" read ul, 
S29, line 13: for 2(4- r) read 2(1 + r) 
S34, line 8: for ug ; read ug; 
S34, line —8: for (0 <x € 1) read (0 < x < 1) 
S54, lines —3, —2: for i read n 
S59, lines — 12, —11: for modulus of the maximum read maximum modulus of the 
S62, line — 3: should read 

since the e's and v's are known and the v's are independent. 
S66, lines —5, —4: for ay, read dys 
S66, line — 1 : should read 

Å = ay Ab + aaf) eu anfa. 
å "be v, Us 


5 
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S69, lines 4 to 7: should read 
Application of Brauer's theorem to this matrix, with a,, taking the possible 
values 1 — 2r, 1 — 2r(1 + h,óx), I — 2r(1 + h;óx), and P, = 2r, shows that 
each of its eigenvalues 2 lies within at least one of the discs 
|4— (1-27 €2r, |A— {1 — 2r(1 + hy dx)}] <2r, A4 — {1 — 2r(1 + h55x))] &2r. 


574, line 12: for mean-value read intermediate value 
S79, line 11 : for on read of 
S98, line —4: delete u, 


E E 
S122, line 6: for J, read Y. 
n=1 


5130, line —6: for M read M, 

S166, line 6: for 70.5 read 70.5? 

S172, solution to exercise 16 : It is unjustifiably assumed that the problem has diagonal 
symmetry in addition to the stated symmetries with respect to Ox and Oy. There 
should therefore be a fourth pivotal value considered, at the point (0, 1), giving a 
fourth equation. The required results can still be deduced though care must be 
taken to ensure that the four equations are consistently ordered. 
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COURSE NOTES Course Notes 


ERRATA IN COURSE UNITS 


Errata in Unit 1 


Page 7, line —3: It is stated that “The MASS of PQ is pAs.” Since any point of the 
string is considered to move only in the transverse direction (page 6, lines —4, — 3) 
the mass of that portion of the string between P and Q must always be the same 
as the mass between points x and x + Ax in the equilibrium position, namely 
pAx. Although the length of this portion of String will certainly increase when in 
motion, its line density will decrease correspondingly as the string stretches. 
p is therefore the equilibrium line density of the string, not the dynamic line density. 


This mistake has repercussions in the ensuing argument. Lines — 2, — / on page? 
are redundant. In line 7 of page 8 the "sec V" should be omitted. In equation U) 
“cos? y” should be “cos”, and in equation (2) "cos* ^ should be "cos? y”. 
Line —9 on page 7 should read 

1 1 i 
se y (1 + tan? på [1 + (8yfoxy]i- 
and “cost yy” should be "cos? ij" in line —7. 


cos? y = 


Erratum in Unit 2 
Page 10, line —3: fort x O read t > 0 


Errata in Unit 3 


Page 7, line —4: for t such t: grad $ read t such that t+ grad à 

Page 8, lines — 18 and — 16: for line — 10 read line —11 

Page 17, line 6: there should be a factor } before the integral sign 

Page 22, line 2: insert (solution on p. 40) underneath this line 

Page 38, line 1: the second integral should be preceded by a minus sign 
Page 42, line 10: for v* read v* 


Errata in Unit 5 


Pages 20, 21, Library Program $CNH321: The version of this program listed is in 
error and has subsequently been corrected and otherwise amended. Those 
interested are invited to obtain a listing of the revised program at their local 
terminal. 

Page 28, line —3: for TË 


= kT, jer read Thi) = KR; 


Errata in Unit 7 


Page I I, line —8: for means read mean 

Page 15, sixth line after the diagram: for v sin 9/AT read v tan O/AT 
Page 19, equation (19): for e *!*! read e "il 

Page 25, line 18: for w read z(x, t) 


DER 1 ay i 1 0% 
: for 3 + 5 read 3 253 
Page 26, line 7: for 2 T ar ead 2 


Page 28, line 7: for ct — x read ct + x 
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